text:
## Call-by-reference decided by caller

I wonder if it makes any sense if the caller of a function decides weither a parameter is passed by reference or value. Often when I look at code (C++ mostly) I ask myself "Is this one changed by the function?" and I have to look up the function declaration. Does anyone know of a language implementing such abstraction?

Edit: I just realized, that the caller can control how an argument is passed by either making a deep copy before passing it (call-by-value) or not (call-by-ref). Only thing is, the function needs to define this argument as ref. Also this does not work well together with manual memory management, because if you create deep copies of stuff on the fly they do not get cleaned up afterwards, because the function receiving the arg does not know it has been copied before.

## Comment viewing options

### ALGOL 60

In ALGOL 60 you can specify whether each individual argument to function is call-by-name or call-by-value, but I believe it is determined in the function definition, not by the caller.

### C# half fixes this

C# doesn't let the caller control whether a parameter is passed by reference or not, but it does require the callsite to indicate it. If you define a method as taking a ref parameter, you must use the ref keyword before the argument when you call it.

### I like this approach. In

I like this approach. In fact recently I heard a lot about C# that sounded very intresting.

### Am I missing something?

Isn't this what you have in C? Parameters are passed by value, if you want to pass by ref, you pass a reference (&s).

### I think I mean sth. slightly

I think I mean sth. slightly different: The caller decides weither a given argument is passed as ref or value, but the function itself treats the parameter as it likes. In C you are bound to the type of parameter defined in the function definition: either a value or a ref (ptr in C).

I just realized that if the caller chooses how to pass an argument, a function may have as well no effect, because its only purpose may be to modify that argument.

### As you say, the definition

As you say, the definition of the function may rely on the parameter being a reference.

You could approach this problem in another way, by introducing an "update with" syntax, for example something like:

foo => myFunction();

This would be syntactic sugar for:

foo = myFunction(foo);

Our functions will always be called by-value, but the caller can decide to replace their old version with this new one immediately. This would also act as a hint to the compiler/interpreter that it can update this value in-place, ie. we can reuse the existing function in a call-by-reference way. Such a hint would also be useful to track argument usage:

myFunction = function(myArg) {
temp = otherFunction(myArg);
return 45 + temp;
}

otherFunction = function(myArg2) {
temp2 = blah(myArg2)
return temp1 * myArg2;
}

Here "myFunction" only uses "myArg" once, so it could do "myArg => otherFunction();" to avoid a copy; this can keep cascading until someone actually needs a copy. An example where a copy is needed is "otherFunction", which uses "myArg2" after defining "temp2", so we have to use two variables in this situation.

This has many similarities to CopyOnWrite; there may be advantages to using both together.

It may be useful to allow name changes without copying too, so we can update values in-place while still giving them more descriptive names. This is especially useful for Hungarian notation. For example, we could change our syntax a little to get:

unsafeEmail = askUser("What is your email address?");
safeEmail => sanitise(unsafeEmail);
sendToDb(safeEmail);

Here we've updated unsafeEmail in-place, but given it a new name "safeEmail". I would recommend that this invalidates the name "unsafeEmail", but there may be arguments for keeping both around (like "$foo =&$bar" in PHP).

### Can be a security issue.

In a language with some kinds of call-by-name, call-by-macro-expansion, call-by-future, or call-by-promise semantics (all of which are uncommon) there is a risk of exposing your caller's environment to examination or even mutation by the callee, so this can become a trust or security issue. The specific problem happens when these dynamic patterns are required to pass a pointer to the caller's environment (so that the argument can be evaluated in that environment) and there is any way for the caller to access that pointer (even via machine code) and see/change what's at the other end of it.

If that is the case, I consider it obligatory for a language design to include a way for the caller to restrict what the callee can see or do. One method of doing this is to provide "restricted environments" for the arguments -- where the argument is reduced to a simple variable name, and the environment passed to the callee with the argument is an environment with only that single variable.

Of course, you'd rather not just explain the technique and then watch as programmers consistently forget to do it when calling routines they have no justification for trusting such as mobile code. Instead, you want to have an explicit "trust" level, probably as part of the "import" or "use" statement, that defines how much trust to extend to a particular callee. This also permits callee functions that *need* access to the caller's environment to work properly, to detect that they do not have it and moves the failure to compile time type-checks rather than runtime.

Of course doing this defeats the purposes behind the evaluation strategy, but there are a whole lot of useful functions that work just fine when called eagerly, and if the callee is one of those it should not matter. If it isn't, then whatever magic it does with your environment simply won't work unless you trust it to mess with your environment.

### Protected from whom?

Most code belongs to the person using it so why do they need protection from themselves? Who can write a C program without passing a pointer to a structure to get something done? In fact, this whole comment seems like a dismissal of all imperative and OOP programming as accessing the data associated with the encapsulated data is an implied pointer from the functions viewpoint.

If you have to worry that somebody with change your data and the language needs to protect you from that person, the person already has access to the programming code and call by reference is the least of your worries.

Of course your solution to this huge problem might just happen to be stateless functional programming. What exactly happens in functional programming when you need the function to act on a 50M data structure? Does a programmer have any rights to the code they write or does it have to pass the compiler writers smell test?

It was nice how you fitted that "compile time type-checks rather than runtime" in there without any justification. It only takes a few words to dismiss dynamic type systems and systems that aren't or can't be compiled monolithically at one time.

So far, all OOP and dynamic languages have bitten the dust. RIP

### Libraries, frameworks, services

Libraries, frameworks, services, plugins, etc. embody a great deal of code not controlled by the person using them. In team development efforts, it is not uncommon that no member of the team understands the full codebase. I wouldn't be surprised if most people control less than 1% of the code they use.

### Calling method determines the code.

If a function is in a library or otherwise doesn't give you access to it's source code, please explain, without changing the function, how to allow the option of call by value or reference? How could this option be available in general through the PL?

You responded to my post and I think I showed that the function would have to be written explicitly to enable both types of calls.

Do you disagree?

### How could this option be

How could this option be available in general through the PL?

The trivial solution would be to compile each function each way that it is used. For functions exposed through an interface, you'd compile them each way. A more sophisticated option is to integrate this with the value model - i.e. a value model that can transparently flip from edit-in-place to copy-on-write (perhaps leveraging uniqueness properties).

The technical feasibility isn't in question. The merit is. I think it would be silly to pass-by-reference where pass-by-value is expected, but the converse might be a little useful.

You responded to my post

I responded to your first sentence. You seem to be depending on a premise (that users own most of their code) that disagrees with my own understandings and experiences.

### I don't think a wrote clearly enough

Let's say I wrote a function to change data that was passed by reference and I return whether it worked or not by using a logical. This function could be arbitrarily complex. If you sent that function a copy of the data instead, the algorithm might work correctly but the result would never get back to the caller. The changed data would just die on the stack.

How would you create an equivalent function that could return the resulting data rather than the logical? What would it mean to the function definition to return 2 totally different types from the same function where that type difference depends on the data? By the way, PHP returns different types for some functions that depend on the data and it is very easy to get errors because of it.

If you can see where my argument is incorrect, please enlighten me!

### There is no general

There is no general equivalence between return values and pass-by-reference values; any such property would depend on the language. But I'm not seeing the topic relevance.

### What exactly happens in

What exactly happens in functional programming when you need the function to act on a 50M data structure?

The usual approach these days is to "stream" the data through our functions, much like we stream text through UNIX pipes. For example I can run grep on my entire hard drive, but it will not wait for all the data before it starts processing.

This isn't completely generic, for example we may have an algorithm which needs access to the entire datastructure before it can do anything, but in those cases we're screwed no matter what language/approach we use.

http://okmij.org/ftp/Streams.html#iteratee

### eh?

Of course your solution to this huge problem might just happen to be stateless functional programming. What exactly happens in functional programming when you need the function to act on a 50M data structure?

[I'll go further than Chris Warburton does in his answer to this]

The equational theory that is the lambda calculus doesn't go all wonky upon reaching a certain length of computation. With that in mind I cannot for the life of me figure out how functional programming at the 50M data structure size is going to be much different than functional programming at any other data structure size.

### Beware the fallacy of the

Beware the fallacy of the beard. Due to various peripheral reasons (performance, process control, partial failure, persistence, memory constraints, limits of human patience), working with large structures is different from working with small ones. Of course, 'large' is relative; these days, 50M isn't too large, but add another order of magnitude or two...

### Indeed, different calling

Indeed, different calling conventions are not equivalent even in a pure, Platonic setting like Lambda Calculus. We usually think of the Church-Rosser property as telling us "all calling conventions give the same result", but we have to remember the caveat "if they give a result at all!". For example:

 infiniteList = 1 : map (1 +) infiniteList

main = print (fst ("Hello world", infiniteList))

If this program halts, it prints "Hello world". This will happen in a non-strict language like Haskell, which won't bother trying to construct the value "infiniteList" because it's not used. A call-by-value language (Scheme, OCaml, etc.) won't give a different result, it will give *no* result; it will get stuck in an infinite loop trying to build the list.

In fact this is the reason I shied away from mentioning calling conventions explicitly; I can't think of a way to make Iteratees for arbitrary, pre-supplied datastructure which is both call-by-value and constant memory. Of course, if we can implement the datastructure too, we can use thunks to generate the large parts on-demand. For example, in Javascript:

 // Naive implementation of the above; will get stuck in an infinite loop. var infiniteList = function() { return [1].concat(infiniteList().map(function(x) { return x + 1; })); }; var main = alert(["Hello world", infiniteList()][0]);

// Implementation using thunks (the second element of each array is a function which generates the rest) var incAll = function(l) { return [1 + l[0], function() { return incAll(l[1]()); }]; }; var infiniteList = function() { return [1, function() { return incAll(infiniteList()); }]; }

If we're getting data from some outside source, these seemingly-artificial thunks can be replaced by the relevant "getNextValue" function for our data source.

### Exactly

Due to various peripheral reasons (performance, process control, partial failure, persistence, memory constraints, limits of human patience),

Exactly, peripheral reasons, none of which are specific to functional programming.

### To counter your earlier

To counter your earlier position, it is sufficient to note that functional programming is different with the larger sizes - different idioms and design patterns, tighter integration with effects for persistence or process-bar feedback, etc.. That this is not specific to functional is a valuable observation for another argument entirely, and perhaps we can study a few paradigms and generalize some requirements for operating on very large data structures.

### Caller determines call by value or reference?

In general you, as a programmer, have access to the function code and you are the one making the call so you could pass by value or reference depending on what you want.

I don't think it makes sense to have the option on the caller side to choose either option. The default in all languages I have used in pass by value and most functions use this automatically. You might override this behavior if you want to change something in your calling space with the function. Nothing wrong with this kind of function, however the code of the function would expect to change your data and then probably return a success/failure code as a return. If you decided to send only a copy of the data then the function behavior would be broken. If the function could detect the kind of variable it was sent then you might be able to make a function that combined both solutions into one function. In general, this kind of complexity wouldn't be worth it.

In the end, it is your code and you can implement it either way without problem and without adding another "feature" to the PL.

### In general you, as a

In general you, as a programmer, have access to the function code and you are the one making the call so you could pass by value or reference depending on what you want.

This statement is incorrect. If I had access to the source in general, that would mean I can always reach the "having-access" state, which I can't (for example, the source may have been lost years ago). In general I *don't* have access, because I can always reach the "no-access" state (eg. by deleting things).

If you mean programmers "usually" have access then I would agree, but hacking library code (without having your patches merged into the trunk) is a bad idea because it puts you in charge of a fork.

### Seems this should be related

Seems this should be somehow related to a problem I was trying to frame a while back, of how to insulate clients from a function's internal notion of side-effects (thereby redressing the global typing problem presented by monads; in a blog post here).

### I say it's all right to require a trust declaration.

I think that if someone wants to be certain that a given module is using only functions that do not have side-effects on anything visible to that module, then it's reasonable for him to say

import-pure "Libraryname1"

Whereas if someone wants to import everything in a given library, even arbitrarily crazy routines that may alter bindings not even related to anything they were passed as arguments, inject arbitrary new bindings into the environment, alter the continuation that the caller will return to when it exits, or do other "crazy" things, it is reasonable for him to be required to say

#I_AM_INSANE //enables absolute trust in imported modules
import-with-insanity "Libraryname2"

At that point, if the compiler detects a call to a routine in Libraryname1 that will have a side effect and modify one of its arguments, it should FAIL because the programmer explicitly told it not to extend such trust to libraryname1. The programmer must find a different way to solve that problem. At the very least the programmer is protected from calling something from libraryname1 which has a visible side effect that he doesn't know about.

Conversely if such a routine were in "Libraryname2" then the compiler wouldn't blink; Libraryname2 has been extended implicit and total trust and may do whatever it wants to the caller and its environment.

### I've been thinking about this because...

In the course of attempting to define a maximally expressive (by Felleisen's notion of expressiveness to be specific) translation-target dialect, it becomes necessary to define ways to express the semantic DESIGN MISTAKES of other languages, because you may be required to translate code from that language, which works only because of those mistakes. At the same time, you'd rather not forego compilation strategies relying on the absence of those mistakes when in fact the program instance in hand does not rely on them.

As others have said, programmers generally use libraries that they do not completely know and/or trust, may not have the source code for, and consider unwise to modify. The problem in this case is to import only those things from that library which are in fact trustworthy. Indeed, you must be able to know (and prove!) to what degree you have limited your code's exposure to the random whimsy that may be expressed in a library originally written in a language that allows the semantic craziness of, for example, macros with arbitrary variable captures.


text:
# How to make a great R reproducible example?

When discussing performance with colleagues, teaching, sending a bug report or searching for guidance on mailing lists and here on SO, a reproducible example is often asked and always helpful.

What are your tips for creating an excellent example? How do you paste data structures from in a text format? What other information should you include?

Are there other tricks in addition to using dput(), dump() or structure()? When should you include library() or require() statements? Which reserved words should one avoid, in addition to c, df, data, etc?

How does one make a great reproducible example?

-
Someone silently voted that this question is off-topic and not in the scope of programming. My response is that this question is very clearly R-specific. I expect the answers to contain references to library, dput and other R-specific code (as in fact they do). Thus I suggest to keep it. –  Andrie May 11 '11 at 17:38
I'm confused about the scope of the question. People seem to have jumped on the interpretation of reproducible example in asking questions on SO or R-help (how to "reproduce the error"). What about reproducible R examples in help pages? In package demos? In tutorials / presentations? –  baptiste Aug 7 '11 at 2:14
@baptiste : The same minus the error. All techniques I explained are used in package help pages, and in tutorials and presentations I give about R –  Joris Meys Oct 4 '11 at 15:11
The data is sometimes the limiting factor, as the structure may be too complex to simulate. To produce public data from private data: stackoverflow.com/a/10458688/742447 in stackoverflow.com/questions/10454973/… –  Etienne Low-Décarie May 11 '12 at 21:41

A minimal reproducible example consists of the following items :

• a minimal dataset, necessary to reproduce the error
• the minimal runnable code necessary to reproduce the error, which can be run on the given dataset.
• the necessary information on the used packages, R version and system it is run on.
• in case of random processes, a seed (set by set.seed()) for reproducibility

Looking at the examples in the help files of the used functions is often helpful. In general, all the code given there fulfills the requirements of a minimal reproducible example : data is provided, minimal code is provided, and everything is runnable.

## Producing a minimal dataset

For most cases, this can be easily done by just providing a vector / dataframe with some values. Or you can use one of the built-in datasets, which are provided with most packages.

Making a vector is easy. Sometimes it is necessary to add some randomness to it, and there are a whole number of functions to make that. sample() can randomize a vector, or give a random vector with only a few values. letters is a useful vector containing the alphabet. This can be used for making factors.

A few examples :

• random values : x <- rnorm(10) for normal distribution, x <- runif(10) for uniform distribution, ...
• a permutation of some values : x <- sample(1:10) for vector 1:10 in random order.
• a random factor : x <- sample(letters[1:4], 20, replace = TRUE)

For matrices, one can use matrix(), eg :

matrix(1:10, ncol = 2)

Making dataframes can be done using data.frame(). One should pay attention to name the entries in the dataframe, and to not make it overly complicated.

An example :

Data <- data.frame(
X = sample(1:10),
Y = sample(c("yes", "no"), 10, replace = TRUE)
)

For some questions, specific formats can be needed. For these, one can use any of the provided as.someType functions : as.factor, as.Date, as.xts, ... These in combination with the vector and/or dataframe tricks.

If you have some data that would be too difficult to construct using these tips, then you can always make a subset of your original data, using eg head(), subset() or the indices. Then use eg. dput() to give us something that can be put in R immediately :

> dput(head(iris,4))
structure(list(Sepal.Length = c(5.1, 4.9, 4.7, 4.6), Sepal.Width = c(3.5,
3, 3.2, 3.1), Petal.Length = c(1.4, 1.4, 1.3, 1.5), Petal.Width = c(0.2,
0.2, 0.2, 0.2), Species = structure(c(1L, 1L, 1L, 1L), .Label = c("setosa",
"versicolor", "virginica"), class = "factor")), .Names = c("Sepal.Length",
"Sepal.Width", "Petal.Length", "Petal.Width", "Species"), row.names = c(NA,
4L), class = "data.frame")

Worst case scenario, you can give a text representation that can be read in using the text parameter of read.table :

zz <- "Sepal.Length Sepal.Width Petal.Length Petal.Width Species
1          5.1         3.5          1.4         0.2  setosa
2          4.9         3.0          1.4         0.2  setosa
3          4.7         3.2          1.3         0.2  setosa
4          4.6         3.1          1.5         0.2  setosa
5          5.0         3.6          1.4         0.2  setosa
6          5.4         3.9          1.7         0.4  setosa"

## Producing minimal code

This should be the easy part, but often isn't. What you should not do, is :

• add all kind of data conversions. Make sure the provided data is already in the correct format (unless that is the problem of course)
• copy-paste a whole function / chunk of code that gives an error. First try to locate which lines exactly result in the error. More often than not you'll find out what the problem is yourself.

What you should do, is :

• add which packages should be used if you use any.
• if you open connections or make files, add some code to close them or delete the files (using unlink())
• if you change options, make sure the code contains a statement to revert them back to the original ones. (eg op <- par(mfrow=c(1,2)) ...some code... par(op) )
• testrun your code in a new, empty R session to make sure the code is runnable. People should be able to just copy-paste your data and your code in the console and get exactly the same as you have.

## Give extra information

In most cases, just the R version and the operating system will suffice. When conflicts arise with packages, giving the output of sessionInfo() can really help. When talking about connections to other applications (be it through ODBC or anything else), one should also provide version numbers for those, and if possible also the necessary information on the setup.

-
How do you use dput if the dataframe is very large and the problem is generated by the middle of the dataframe? Is there a way to use dput to reproduce the mid-section of data, say rows 60 through 70? –  BgnR Apr 10 '14 at 13:12
@BgnR You can extract part of the data frame using indices, eg: tmp <- mydf[50:70,] followed by dput(mydf). If the data frame is really big, try isolating the problem and just submit the few lines that cause the problem. –  Joris Meys Apr 11 '14 at 11:27
@JorisMeys: Is there a way to tell head or dput to limit data to level N recursively? I'm trying to come up with reproducible example and my data is a list of data frames. So, dput(head(myDataObj)) seems not to be enough, as it generates an output file of 14MB size. –  Aleksandr Blekh Aug 4 '14 at 15:51
@JorisMeys: Just FYI - posted question in the comment above as a separate question: stackoverflow.com/questions/25127026/…. –  Aleksandr Blekh Aug 4 '14 at 20:32

# How to write a reproducible example.

You are most likely to get good help with your R problem if you provide a reproducible example. A reproducible example allows someone else to recreate your problem by just copying and pasting R code.

There are four things you need to include to make your example reproducible: required packages, data, code, and a description of your R environment.

• Packages should be loaded at the top of the script, so it's easy to see which ones the example needs.

• The easiest way to include data in an email or Stack Overflow question is to use dput() to generate the R code to recreate it. For example, to recreate the mtcars dataset in R, I'd perform the following steps:

1. Run dput(mtcars) in R
2. Copy the output
3. In my reproducible script, type mtcars <- then paste.
• Spend a little bit of time ensuring that your code is easy for others to read:

• make sure you've used spaces and your variable names are concise, but informative

• do your best to remove everything that is not related to the problem.
The shorter your code is, the easier it is to understand.

• Include the output of sessionInfo() in a comment in your code. This summarises your R environment and makes it easy to check if you're using an out-of-date package.

You can check you have actually made a reproducible example by starting up a fresh R session and pasting your script in.

Before putting all of your code in an email, consider putting it on http://gist.github.com/. It will give your code nice syntax highlighting, and you don't have to worry about anything getting mangled by the email system.

-

Personally, I prefer "one" liners. Something along the lines:

my.df <- data.frame(col1 = sample(c(1,2), 10, replace = TRUE),
col2 = as.factor(sample(10)), col3 = letters[1:10],
col4 = sample(c(TRUE, FALSE), 10, replace = TRUE))
my.list <- list(list1 = my.df, list2 = my.df[3], list3 = letters)

The data structure should mimic the idea of writer's problem, and not the exact verbatim structure. I really appreciate it when variables don't overwrite my own variables or god forbid, functions (like df).

Alternatively, one could cut a few corners and point to a pre-existing data set, something like:

library(vegan)
data(varespec)
ord <- metaMDS(varespec)

Don't forget to mention any special packages you might be using.

If you're trying to demonstrate something on larger objects, you can try

my.df2 <- data.frame(a = sample(10e6), b = sample(letters, 10e6, replace = TRUE))

If you're working with spatial data via the raster package, you can generate some random data. A lot of examples can be found in the package vignette, but here's a small nugget.

library(raster)
r1 <- r2 <- r3 <- raster(nrow=10, ncol=10)
values(r1) <- runif(ncell(r1))
values(r2) <- runif(ncell(r2))
values(r3) <- runif(ncell(r3))
s <- stack(r1, r2, r3)

If you're in need of some spatial object as implemented in sp, you can get some data sets via external files (like ESRI shape file) in "spatial" packages (see Spatial view in Task Views).

library(rgdal)
ogrDrivers()
dsn <- system.file("vectors", package = "rgdal")[1]
ogrListLayers(dsn)
ogrInfo(dsn=dsn, layer="cities")

-
+1 for the built-in datasets. –  Joris Meys May 11 '11 at 11:46

Inspired by this very post, I now use a handy function
reproduce(<mydata>) when I need to post to StackOverflow.

## QUICK INSTRUCTIONS

If myData is the name of your object to reproduce, run the following in R:

install.packages("devtools")
library(devtools)
source_url("https://raw.github.com/rsaporta/pubR/gitbranch/reproduce.R")

reproduce(myData)

## Details:

This function is an intelligent wrapper to dput and does the following:

• automatically samples a large data set (based on size and class. Sample size can be adjusted)
• creates a dput output
• allows you to specify which columns to export
• appends to the front of it objName <- ... so that it can be easily copy+pasted, but...
• If working on a mac, the output is automagically copied to the clipboard, so that you can simply run it and then paste to your question.

## Example:

# sample data
DF <- data.frame(id=rep(LETTERS, each=4)[1:100], replicate(100, sample(1001, 100)), Class=sample(c("Yes", "No"), 100, TRUE))

DF is about 100 x 102. I want to sample 10 rows, and a few specific columns

reproduce(DF, cols=c("id", "X1", "X73", "Class"))  # I could also specify the column number.

### Gives the following output:

This is what the sample looks like:

id  X1 X73 Class
1    A 266 960   Yes
2    A 373 315    No            Notice the selection split
3    A 573 208    No           (which can be turned off)
4    A 907 850   Yes
5    B 202  46   Yes
6    B 895 969   Yes   <~~~ 70 % of selection is from the top rows
7    B 940 928    No
98   Y 371 171   Yes
99   Y 733 364   Yes   <~~~ 30 % of selection is from the bottom rows.
100  Y 546 641    No

==X==============================================================X==
Copy+Paste this part. (If on a Mac, it is already copied!)
==X==============================================================X==

DF <- structure(list(id = structure(c(1L, 1L, 1L, 1L, 2L, 2L, 2L, 25L, 25L, 25L), .Label = c("A", "B", "C", "D", "E", "F", "G", "H", "I", "J", "K", "L", "M", "N", "O", "P", "Q", "R", "S", "T", "U", "V", "W", "X", "Y"), class = "factor"), X1 = c(266L, 373L, 573L, 907L, 202L, 895L, 940L, 371L, 733L, 546L), X73 = c(960L, 315L, 208L, 850L, 46L, 969L, 928L, 171L, 364L, 641L), Class = structure(c(2L, 1L, 1L, 2L, 2L, 2L, 1L, 2L, 2L, 1L), .Label = c("No", "Yes"), class = "factor")), .Names = c("id", "X1", "X73", "Class"), class = "data.frame", row.names = c(1L, 2L, 3L, 4L, 5L, 6L, 7L, 98L, 99L, 100L))

==X==============================================================X==

Notice also that the entirety of the output is in a nice single, long line, not a tall paragraph of chopped up lines. This makes it easier to read on SO questions posts and also easier to copy+paste.

# ______________________________________

## Update Oct 2013:

You can now specify how many lines of text output will take up (ie, what you will paste into StackOverflow). Use the lines.out=n argument for this. Example:

reproduce(DF, cols=c(1:3, 17, 23), lines.out=7) yields:

    ==X==============================================================X==
Copy+Paste this part. (If on a Mac, it is already copied!)
==X==============================================================X==

DF <- structure(list(id = structure(c(1L, 1L, 1L, 1L, 2L, 2L, 2L, 25L,25L, 25L), .Label
= c("A", "B", "C", "D", "E", "F", "G", "H","I", "J", "K", "L", "M", "N", "O", "P", "Q", "R", "S", "T", "U","V", "W", "X", "Y"), class = "factor"),
X1 = c(809L, 81L, 862L,747L, 224L, 721L, 310L, 53L, 853L, 642L),
X2 = c(926L, 409L,825L, 702L, 803L, 63L, 319L, 941L, 598L, 830L),
X16 = c(447L,164L, 8L, 775L, 471L, 196L, 30L, 420L, 47L, 327L),
X22 = c(335L,164L, 503L, 407L, 662L, 139L, 111L, 721L, 340L, 178L)), .Names = c("id","X1",
"X2", "X16", "X22"), class = "data.frame", row.names = c(1L,2L, 3L, 4L, 5L, 6L, 7L, 98L, 99L, 100L))

==X==============================================================X==

-
+1, I like the function, but for me your closing paragraph is a bit subjective. As someone who looks at a lot of R code, I would rather see a tall paragraph where I can get a sense of what I'll be creating, and that's much easier to do if it is all there on my screen, and not in a format where I need to scroll horizontally. –  Ananda Mahto Oct 24 '13 at 4:11
@AnandaMahto, you have a great point. I've added a lines.out argument which allows you to break up the output. The function takes a best guess at logical breaks, but this is purely regex based. Convceivably, we could grep the output (before it's returned) for names of the data and break there. Feel free to look at the source and edit if you'd like! –  Ricardo Saporta Oct 24 '13 at 5:03
@AnandaMahto, the function now uses the names of the value (if available) to break out the lines. –  Ricardo Saporta Oct 24 '13 at 6:03

Here is a good guide:

http://www.r-bloggers.com/three-tips-for-posting-good-questions-to-r-help-and-stack-overflow/

But the most important is: Just make sure that you make a small piece of code that we can run to see what the problem is. A usefull function for this is dput(), but if you have very large data you might want to make a small sample dataset or only use the first 10 lines or so.

EDIT:

Also make sure that you identified where the problem is yourself. The example should not be an entire R script with "On line 200 there is an error". If you use the debugging tools in R (I love browser()) and google you should be able to really identify where the problem is and reproduce a trivial example in which the same thing goes wrong.

-

The R-help mailing list has a posting guide which covers both asking and answering questions, including an example of generating data:

Examples: Sometimes it helps to provide a small example that someone can actually run. For example:

If I have a matrix x as follows:

  > x <- matrix(1:8, nrow=4, ncol=2,
dimnames=list(c("A","B","C","D"), c("x","y"))
> x
x y
A 1 5
B 2 6
C 3 7
D 4 8
>

how can I turn it into a dataframe with 8 rows, and three columns named 'row', 'col', and 'value', which have the dimension names as the values of 'row' and 'col', like this:

  > x.df
row col value
1    A   x      1

...
(To which the answer might be:

  > x.df <- reshape(data.frame(row=rownames(x), x), direction="long",
varying=list(colnames(x)), times=colnames(x),
v.names="value", timevar="col", idvar="row")

)

The word small is especially important. You should be aiming for a minimal reproducible example, which means that the data and the code should be as simple as possible to explain the problem.

EDIT: Pretty code is easier to read than ugly code. Use a style guide.

-

Since R.2.14 (I guess) you can feed your data text representation directly to read.table:

df <- read.table(header=T, text="Sepal.Length Sepal.Width Petal.Length Petal.Width Species
1          5.1         3.5          1.4         0.2  setosa
2          4.9         3.0          1.4         0.2  setosa
3          4.7         3.2          1.3         0.2  setosa
4          4.6         3.1          1.5         0.2  setosa
5          5.0         3.6          1.4         0.2  setosa
6          5.4         3.9          1.7         0.4  setosa
")

-
You can also use read.table("clipboard", header=TRUE). –  sebastian-c Sep 17 '12 at 1:09
@sebastian-c how is that good for making reproducible example?? :) –  TMS May 24 '14 at 8:43
@TMS Giving it serious thought, if the asker has supplied the data and the problem is small (but might have a few solutions), then it might be faster and you can still follow all the steps. –  sebastian-c May 26 '14 at 17:52

Sometimes the problem really isn't reproducible with a smaller piece of data, no matter how hard you try, and doesn't happen with synthetic data (although it's useful to show how you produced synthetic data sets that did not reproduce the problem, because it rules out some hypotheses).

• Posting the data to the web somewhere and providing a URL may be necessary.
• If the data can't be released to the public at large but could be shared at all, then you may be able to offer to e-mail it to interested parties (although this will cut down the number of people who will bother to work on it).
• I haven't actually seen this done, because people who can't release their data are sensitive about releasing it any form, but it would seem plausible that in some cases one could still post data if it were sufficiently anonymized/scrambled/corrupted slightly in some way.

If you can't do either of these then you probably need to hire a consultant to solve your problem ...

edit: Two useful SO questions for anonymization/scrambling:

-
For producing synthetic data sets, the answers to this question give useful examples, including applications of fitdistr and fitdistrplus. –  Iterator Oct 23 '11 at 13:47

The answers so far are obviously great for the reproducibility part. This is merely to clarify that a reproducible example cannot and should not be the sole component of a question. Don't forget to explain what you want it to look like and the contours of your problem, not just how you have attempted to get there so far. Code is not enough; you need words also.

Here's a reproducible example of what to avoid doing (drawn from a real example, names changed to protect the innocent):

The following is sample data and part of function I have trouble with.

code
code
code
code
code (40 or so lines of it)

How can I achieve this ?

-

To quickly create a dput of your data you can just copy (a piece of) the data to your clipboard and run the following in R:

for data in Excel:

dput(read.table("clipboard",sep="\t",header=TRUE))

for data in a txt file:

dput(read.table("clipboard",sep="",header=TRUE))

You can change the sep in the latter if necessary. This will only work if your data is in the clipboard of course.

-

Often you need some data for an example, however, you don't want to post your exact data. To use some existing data.frame in established library, use data command to import it.

e.g.,

data(mtcars)

and then do the problem

names(mtcars)
your problem demostrated on the mtcars data set

-
Many built-in data sets (like popular mtcars and iris datasets) don't actually need the data call to be used. –  Gregor Jul 17 '13 at 19:17

If you have large dataset which cannot be easily put to the script using dput(), post your data to pastebin and load them using read.table:

d <- read.table("http://pastebin.com/raw.php?i=m1ZJuKLH")

Inspired by @Henrik.

-

I have a very easy and efficient way to make a R example that has not been mentioned above. You can define your structure firstly.For example,

mydata <- data.frame(a=character(0), b=numeric(0),  c=numeric(0), d=numeric(0))

>fix(mydata)

Then you can input your data manually.This is efficient for smaller examples rather than big ones.

-
...then dput(mydata) –  GSee Mar 8 '14 at 16:46

Reproducible code is key to get help. However, there are many users that might be skeptical of pasting even a chunk of their data. For instance, they could be working with sensitive data or on an original data collected to use in a research paper. For any reason, I thought it would be nice to have a handy function for "deforming" my data before pasting it publicly. The anonymize function from the package SciencesPo is very silly, but for me it works nicely with dput function.

install.packages("SciencesPo")

dt <- data.frame(
Z = sample(LETTERS,10),
X = sample(1:10),
Y = sample(c("yes", "no"), 10, replace = TRUE)
)

> dt
Z  X   Y
1  D  8  no
2  T  1 yes
3  J  7  no
4  K  6  no
5  U  2  no
6  A 10 yes
7  Y  5  no
8  M  9 yes
9  X  4 yes
10 Z  3  no

Then I anonymize it:

> anonymize(dt)
Z    X  Y
1   b2  2.5 c1
2   b6 -4.5 c2
3   b3  1.5 c1
4   b4  0.5 c1
5   b7 -3.5 c1
6   b1  4.5 c2
7   b9 -0.5 c1
8   b5  3.5 c2
9   b8 -1.5 c2
10 b10 -2.5 c1

One may also want to sample few variables instead of the whole data before apply anonymization and dput command.

    # sample two variables without replacement
> anonymize(sample.df(dt,5,vars=c("Y","X")))
Y    X
1 a1 -0.4
2 a1  0.6
3 a2 -2.4
4 a1 -1.4
5 a2  3.6

-

I wonder if an http://www.r-fiddle.org/ link could be a very neat way of sharing a problem. It receives a unique ID like and one could even think about embedding it in SO.

http://www.r-fiddle.org/#/help

-

If you have one or more factor variable(s) in your data that you want to make reproducible with dput(head(mydata)), consider adding droplevels to it, so that levels of factors that are not present in the minimized data set are not included in your dput output, in order to make the example minimal:

dput(droplevels(head(mydata)))

-

## protected by GSeeMay 31 '14 at 2:51

Thank you for your interest in this question. Because it has attracted low-quality answers, posting an answer now requires 10 reputation on this site.

format: plain


text:
# Ch. 2 Statistical Hypothesis Testing

## 2.1 Hypotheses and Test-Statistics

Assume an independently and identically distributed (i.i.d.) random sample $$X_1,\dots,X_n$$, where the distribution of $$X_i$$, $$i=1,\dots,n$$, depends on some unknown parameter $$\theta\in\Omega$$, and where $$\Omega$$ is some parameter space.

General Testing Problem:

$H_0:\theta\in\Omega_0$ against $H_1:\theta\in\Omega_1$

$$H_0$$ is the null hypothesis, while $$H_1$$ is the alternative. $$\Omega_0\subset\Omega$$ and $$\Omega_1\subset\Omega$$ are used to denote the possible values of $$\theta$$ under $$H_0$$ and $$H_1$$. Necessarily, $$\Omega_0\cap\Omega_1=\emptyset$$.

For a large number of tests we have $$\Omega\subseteq\mathbb{R}$$ and the respective null hypothesis states that $$\theta$$ has a specific value $$\theta_0\in\mathbb{R}$$, i.e., $$\Omega_0=\{\theta_0\}$$ and $$H_0:\theta=\theta_0$$. Depending on the alternative one then often distinguishes between one-sided ($$\Omega_1=(\theta_0,\infty)$$ or $$\Omega_1=(-\infty,\theta_0)$$) and two-sided tests ($$\Omega_1=\{\theta\in\mathbb{R}|\theta\neq \theta_0\}$$).

The data $$X_1,\dots,X_n$$ is used in order to decide whether to accept or to reject $$H_0$$.

Test Statistic: Every statistical hypothesis test relies on a corresponding test statistic $T=T(X_1,\dots,X_n).$ Any test statistic is a real valued random variable, and for given data the resulting observed value $$T_{obs}$$ is used to decide between $$H_0$$ and $$H_1$$. Generally, the distribution of $$T$$ under $$H_0$$ is analyzed in order to define a rejection region $$C$$:

• $$T_{obs}\not\in C$$ $$\Rightarrow$$ $$H_0$$ is not rejected
• $$T_{obs}\in C$$ $$\Rightarrow$$ $$H_0$$ is rejected

For one-sided tests $$C$$ is typically of the form $$(-\infty,c_0]$$ or $$[c_1,\infty)$$. For two-sided tests $$C$$ typically takes the form of $$(-\infty,c_0]\cup [c_1,\infty)$$. The limits $$c_0$$ and $$c_1$$ of the respective intervals are called critical values, and are obtained from quantiles of the null distribution, i.e., the distribution of $$T$$ under $$H_0$$.

Decision Errors:

Decision Errors Verbal Definition Probability of a Type I/II Error
Type I error $$H_0$$ is rejected even though $$H_0$$ is true. $$P(T\in C |\,H_0\text{ true})$$
Type II error The test fails to reject a false $$H_0$$. $$P(T\not\in C |\,H_1\text{ true})$$

## 2.2 Significance Level, Size and p-Values

Significance Level: In a statistical significance test, the probability of a type I error is controlled by the significance level $$\alpha$$ (e.g., $$\alpha=5\%$$).

$P\left(\text{Type I error}\right)=\sup_{\theta\in\Omega_0} P(T\in C|\theta\in\Omega_0)\leq \alpha$

Size: The size of a statistical test is defined as $\sup_{\theta\in\Omega_0} P(T\in C|\theta\in\Omega_0).$

That is, the preselected significance level $$\alpha$$ is an upper bound for the size, which may not be attained (i.e., size $$<\alpha$$) if, for instance, the relevant probability function is discrete.

Practically important significance levels:

• $$\alpha=0.05$$: It is common to say that a test result is “significant” if a hypothesis test of level $$\alpha=0.05$$ rejects $$H_0$$.
• $$\alpha=0.01$$: It is common to say that a test result is “strongly significant” if a hypothesis test of level $$\alpha=0.01$$ rejects $$H_0$$.

p-Value: The p-value is the probability of obtaining a test statistic at least as “extreme” as the one that was actually observed, assuming that the null hypothesis is true.

• For one-sided tests:
• $$P(T\geq T_{\text{obs}}|H_0\text{ true})$$ or
• $$P(T\leq T_{\text{obs}}|H_0\text{ true})$$
• For two-sided tests:
• $$2\min\{P(T\leq T_{\text{obs}}|H_0\text{ true}),\,P(T\geq T_{\text{obs}}|H_0\text{ true})\}$$

Remarks:

• The p-value is random as it depends on the observed data. That is, different random samples will lead to different p-values.
• For given data, having determined the p-value of a test we also know the test decisions for all possible levels $$\alpha$$:

• $$\alpha > \text{p-value} \Rightarrow H_0 \text{ is rejected}$$
• $$\alpha < \text{p-value} \Rightarrow H_0 \text{ cannot be rejected}$$

Example: Let $$X_i\sim N(\mu,\sigma^2)$$ independently for all $$i=1,\dots,5=n$$. Observed realizations from this i.i.d. random sample: $$X_1=19.20$$, $$X_2=17.40$$, $$X_3=18.50$$, $$X_4=16.50$$, $$X_5=18.90$$. That is, the empirical mean is given by $$\bar X =18.1$$.

Testing problem: $$H_0:\mu=\mu_0$$ against $$H_1:\mu\ne\mu_0$$ (i.e., a two-sided test), where $$\mu_0=17$$.

Since the variance is unknown, we use the sample standard deviation, $$s$$, which then leads to the t-test for testing $$H_0$$. Test statistic of the t-test: $T=\frac{\sqrt{n}(\bar X-\mu_0)}{s},$ where $$s^2=\frac{1}{n-1}\sum_{i=1}^n (X_i-\bar X)^2$$ is the unbiased estimator of $$\sigma^2$$. $T_{obs}=\frac{\sqrt{5}(18.1-17)}{1.125}=2.187$ $\Rightarrow \hbox{p-value}=2\min\{P(T_{n-1}\leq 2.187),\, P(T_{n-1}\geq 2.187)\}=0.094$

The above computations in R

library("magrittr", quietly = TRUE)# for using the pipe-operator: %>%

X           <- c(19.20, 17.40, 18.50, 16.50, 18.90)
mu_0        <- 17        # hypothetical mean
n           <- length(X) # sample size
X_mean      <- mean(X)   # empirical mean
X_sd        <- sd(X)     # empirical sd
# t-test statistic
t_test_stat <- sqrt(n)*(X_mean - mu_0)/X_sd

# p-value for two-sided test
c(pt(q = t_test_stat, df = n-1, lower.tail = TRUE),
pt(q = t_test_stat, df = n-1, lower.tail = FALSE)) %>%
min * 2 -> p_value

p_value %>% round(., digits = 3)
## [1] 0.094

Of course, there is also a t.test() function in R:

t.test(X, mu = mu_0, alternative = "two.sided")
##
##  One Sample t-test
##
## data:  X
## t = 2.1869, df = 4, p-value = 0.09402
## alternative hypothesis: true mean is not equal to 17
## 95 percent confidence interval:
##  16.70347 19.49653
## sample estimates:
## mean of x
##      18.1

## 2.3 The Power Function

For every possible value $$\theta\in\Omega_0\cup\Omega_1$$, all sample sizes $$n$$ and each significance level $$\alpha$$ the corresponding value of the power function $$\beta$$ is defined by the following probability: $\beta_{n,\alpha}(\theta):=P(\text{reject } H_0\; |\; \theta\in\Omega_0\cup\Omega_1)$

Obviously, $$\beta_{n,\alpha}(\theta)\leq \alpha$$ for all $$\theta\in\Omega_0$$. Furthermore, for any $$\theta\in\Omega_1$$, $$1-\beta_{n,\alpha}(\theta)$$ is the probability of committing a type II error.

The power function is an important tool for accessing the quality of a test and for comparing different test procedures.

Conservative Test: If possible, a test is constructed in such a way that size equals level, i.e., $$\beta_{n,\alpha}(\theta)=\alpha$$ for some $$\theta\in\Omega_0$$. In some cases, however, as for discrete test statistics or complex, composite null hypothesis, it is not possible to reach the level, and $$\sup_{\theta\in\Omega_0}\beta_{n,\alpha}(\theta)<\alpha$$. In this case the test is called conservative.

Unbiased Test: A significance test of level $$\alpha>0$$ is called unbiased if $$\beta_{n,\alpha}(\theta)\ge\alpha$$ for all $$\theta\in\Omega_1$$.

Consistent Test: A significance test of level $$\alpha>0$$ is called consistent if $\lim_{n\rightarrow \infty} \beta_{n,\alpha}(\theta) =1$ for all $$\theta\in \Omega_1$$.

Most Powerful Test: When choosing between different testing procedures for the same testing problem, one will usually prefer the most powerful test. Consider a fixed sample size $$n$$. For a specified $$\theta\in\Omega_1$$, a test with power function $$\beta_{n,\alpha}(\theta)$$ is said to be most powerful for $$\theta$$ if for any alternative test with power function $$\beta^*_{n,\alpha}(\theta)$$, $\beta_{n,\alpha}(\theta)\ge \beta^*_{n,\alpha}(\theta)$ holds for all levels $$\alpha>0$$.

Uniformly Most Powerful: A test with power function $$\beta_{n,\alpha}(\theta)$$ is said to be uniformly most powerful against the set of alternatives $$\Omega_1$$ if for any alternative test with power function $$\beta^*_{n,\alpha}(\theta)$$, $\beta_{n,\alpha}(\theta)\ge \beta^*_{n,\alpha}(\theta)\quad \text{holds for all }\theta\in\Omega_1, \alpha>0$ Unfortunately, uniformly most powerful tests only exist for very special testing problems.

Example: Let $$X_1,\dots,X_n$$ be an i.i.d. random sample. Assume that $$n=9$$, and that $$X_i\sim N(\mu,0.18^2)$$. Hence, in this simple example only the mean $$\mu=E(X)$$ is unknown, while the standard deviation has the known value $$\sigma=0.18$$.

Testing problem: $$H_0:\mu=\mu_0$$ against $$H_1:\mu\neq \mu_0$$ for $$\mu_0=18.3$$ (i.e., a two-sided test).

Since the variance is known, a test may rely on the Gauss (or Z) test statistic: $Z=\frac{\sqrt{n} (\bar X -\mu_0)}{\sigma} =\frac{3 (\bar X -18.3)}{0.18}$

Under $$H_0$$ we have $$Z\sim N(0,1)$$, and for the significance level $$\alpha=0.05$$ the null hypothesis is rejected if $|Z|\geq z_{1-\alpha/2}=1.96,$ where $$z_{1-\alpha/2}$$ denotes the $$(1-\alpha/2)$$-quantile of the standard normal distribution. Note that the size of this test equals its level $$\alpha=0.05$$.

For determining the rejection region of a test it suffices to determine the distribution of the test statistic under $$H_0$$. But in order to calculate the power function one needs to quantify the distribution of the test statistic for all possible values $$\theta\in\Omega$$. For many important problems this is a formidable task. For the Gauss test, however, it is quite easy. Note that for any (true) mean value $$\mu\in\mathbb{R}$$ the corresponding distribution of $$Z\equiv Z_\mu=\sqrt{n(\bar X-\mu_0)}/\sigma$$ is $Z_\mu=\frac{\sqrt{n} (\mu -\mu_0)}{\sigma}+\frac{\sqrt{n} (\bar X -\mu)}{\sigma} \sim N\left(\frac{\sqrt{n} (\mu -\mu_0)}{\sigma}, 1\right)$ This implies that \begin{align*} \beta_{n,\alpha}(\mu) & = P\left(|Z_\mu|>z_{1-\alpha/2}\right)\\ & = 1-\Phi\left(z_{1-\alpha/2}-\frac{\sqrt{n} (\mu -\mu_0)}{\sigma}\right) + \Phi\left(-z_{1-\alpha/2}-\frac{\sqrt{n} (\mu -\mu_0)}{\sigma}\right), \end{align*} where $$\Phi$$ denotes the distribution function of the standard normal distribution.

Implementing the power function of the two-sided Z-test in R:

# The power function
beta_Ztest_TwoSided <- function(n, alpha, sigma, mu_0, mu){
# (1-alpha/2)-quantile of N(0,1):
z_upper        <- qnorm(p = 1-alpha/2)
# location shift under H_1:
location_shift <- sqrt(n) * (mu - mu_0)/sigma
# compute power
power          <- 1 - pnorm( z_upper - location_shift) +
pnorm(-z_upper - location_shift)
return(power)
}

# Apply the function
n     <-  9
sigma <-  0.18
mu_0  <- 18.3
##
c(beta_Ztest_TwoSided(n = n, alpha = 0.05, sigma = sigma, mu_0 = mu_0, mu=18.35),
beta_Ztest_TwoSided(n = n, alpha = 0.05, sigma = sigma, mu_0 = mu_0, mu=18.50),
beta_Ztest_TwoSided(n = n, alpha = 0.01, sigma = sigma, mu_0 = mu_0, mu=18.50)) %>%
round(., digits = 3)
## [1] 0.133 0.915 0.776

Plotting the graph of the power function

suppressPackageStartupMessages(
library("tidyverse")
)
# Vectorize the function with respect to mu_0:
beta_Ztest_TwoSided <- Vectorize(FUN = beta_Ztest_TwoSided,
vectorize.args = "mu_0")

mu_0_vec <- seq(from = 17.75, to = 18.25, len = 50)

beta_vec <- beta_Ztest_TwoSided(n     =   10,
alpha = 0.05,
sigma = 0.18,
mu    =  18,
mu_0  = mu_0_vec)

beta_df <- data.frame("mu_0"  = mu_0_vec,
"Beta"  = beta_vec)

ggplot(data = beta_df, aes(x=mu_0, y=Beta)) +
geom_line() +
geom_hline(yintercept = 0.05, lty=2) +
geom_text(aes(x=17.77, y=0.07, label='alpha==0.05'), parse=TRUE, size=5) +
labs(title = expression(
paste("Powerfunction of the two-sided Z-Test (n=10 and ",alpha==0.05,")")),
x = expression(paste(mu[0])),
y = expression(paste(beta)), size=8)    +
theme_bw() +
theme(axis.text  = element_text(size=12),
axis.title = element_text(size=14))

This example illustrates the power function of a sensible test, since:

• Under $$H_0:\mu=\mu_0$$ we have $$\beta_{n,\alpha}(\mu_0)=\alpha$$.
• The test is unbiased, since $$\beta_{n,\alpha}(\mu)\geq\alpha$$ for any $$\mu\neq\mu_0$$.
• The test is consistent, since $$\lim_{n\rightarrow\infty} \beta_{n,\alpha}(\mu)=1$$ for every fixed $$\mu\neq \mu_0$$.
• For fixed sample size $$n$$, $$\beta_{n,\alpha}(\mu)$$ increases as the distance $$|\mu-\mu_0|$$ increases.
• If $$|\mu-\mu_0|>|\mu^*-\mu_0|$$ then $$\beta_{n,\alpha}(\mu)>\beta_{n,\alpha}(\mu^*)$$.
• $$\beta_{n,\alpha}(\mu)$$ decreases as the significance level $$\alpha$$ of the test decreases. I.e., if $$\alpha>\alpha^*$$ then $$\beta_{n,\alpha}(\mu)>\beta_{n,\alpha^*}(\mu)$$.

Assuming that the basic assumptions (i.e., normality and known variance) are true, the above Gauss-test is the most prominent example of a uniformly most powerful test. Under its (restrictive) assumptions, no other possible test can achieve a larger value of $$\beta_{n,\alpha}(\mu)$$ for any possible value of $$\mu$$.

## 2.4 Asymptotic Null Distributions

Generally, the underlying distributions are unknown. In this case it is usually not possible to compute the power function of a test for fixed $$n$$. (Exceptions are so called “distribution-free” tests in nonparametric statistics.) The only way out of this difficulty is to rely on large sample asymptotics and corresponding asymptotic distributions, which allow to approximate the power function and to study the asymptotic efficiency of a test. The finite sample behavior of a test for different sample sizes $$n$$ is then evaluated by means of simulation studies.

For a real-valued parameter $$\theta$$ most tests of $$H_0:\theta=\theta_0$$ rely on estimators $$\hat\theta$$ of $$\theta$$. Under suitable regularity conditions on the underlying distribution, central limit theorems usually imply that $\sqrt{n}(\hat\theta - \theta)\rightarrow_D N(0,v^2)\quad\text{as}\quad n\rightarrow\infty,$ where $$v^2$$ is the asymptotic variance of the estimator.

Often a consistent estimator $$\hat v^2$$ of $$v^2$$ can be determined from the data. For large $$n$$ we then approximately have $\frac{\sqrt{n}(\hat\theta - \theta)}{ v}\overset{a}{\sim} N(0,1).$ For a given $$\alpha$$, a one-sided test of $$H_0:\theta=\theta_0$$ against $$H_1:\theta>\theta_0$$ then rejects $$H_0$$ if $Z=\frac{\sqrt{n} (\hat\theta -\theta_0)}{v}>z_{1-\alpha}.$ The corresponding asymptotic approximation (valid for sufficiently large $$n$$) of the true power function is then given by $\beta_{n,\alpha}(\theta) = 1-\Phi\left(z_{1-\alpha}-\frac{\sqrt{n} (\theta -\theta_0)}{v}\right)$

Note that in practice the (unknown) true value $$v^2$$ is generally replaced by an estimator $$\hat v^2$$ determined from the data. As long as $$\hat v^2$$ is a consistent estimator of $$v^2$$ this leads to the same asymptotic power function. The resulting test is asymptotically unbiased and consistent.

Usually there are many different possible estimators for a parameter $$\theta$$. Consider an alternative estimator $$\tilde\theta$$ of $$\theta$$ satisfying $\sqrt{n}(\tilde\theta - \theta)\rightarrow_D N(0,\tilde v^2) \quad\text{as}\quad n\rightarrow\infty.$ If the asymptotic variance $$v^2$$ of the estimator $$\hat\theta$$ is smaller than the asymptotic variance $$\tilde v^2$$ of $$\tilde\theta$$, i.e., $$v^2<\tilde v^2$$, then $$\hat\theta$$ is a more efficient estimator of $$\theta$$. Then necessarily the test based on $$\hat\theta$$ is more powerful than the test based on $$\tilde\theta$$, since asymptotically for all $$\theta>\theta_0$$ \begin{align*} \tilde\beta_{n,\alpha}(\theta) &= 1-\Phi\left(z_{1-\alpha}-\frac{\sqrt{n} (\theta -\theta_0)}{\tilde v}\right)\\ & < 1-\Phi\left(z_{1-\alpha}-\frac{\sqrt{n} (\theta -\theta_0)}{v}\right)=\beta_{n,\alpha}(\theta) \end{align*}

Example: Let $$X_1,\dots,X_n$$ be an iid random sample. Consider testing $$H_0:\mu=\mu_0$$ against $$H_1:\mu>\mu_0$$, where $$\mu:=E(X_i)$$. For a given level $$\alpha$$ the t-test then rejects $$H_0$$ if $T=\frac{\sqrt{n}(\bar X-\mu_0)}{S}>t_{n-1;1-\alpha},$ where $$t_{n-1;1-\alpha}$$ is the $$1-\alpha$$ quantile of a t-distributions with $$n-1$$-degrees of freedom. This is an exact test if the distribution of $$X_i$$ is normal. In the general case, the justification of the t-test is based on asymptotic arguments. Under some regularity conditions the central limit theorem implies that $\sqrt{n}(\bar X - \mu)\rightarrow_D N(0,\sigma^2)\quad\text{as}\quad n\rightarrow\infty$ with $$\sigma^2=Var(X_i)$$. Moreover, $$S^2$$ is a consistent estimator of $$\sigma^2$$ and $$t_{n-1;1-\alpha}\rightarrow z_{1-\alpha}$$ as $$n\rightarrow \infty$$. Thus even if the distribution of $$X_i$$ is non-normal, for sufficiently large $$n$$, $$T=\frac{\sqrt{n}(\bar X-\mu_0)}{S}$$ is approximately $$N(0,1)$$-distributed and the asymptotic power function of the t-test is given by $\beta_{n,\alpha}(\theta) = 1-\Phi\left(z_{1-\alpha}-\frac{\sqrt{n} (\mu -\mu_0)}{\sigma}\right).$

## 2.5 Multiple Comparisons

In statistics, the multiple comparisons, multiplicity or multiple testing problem occurs when one considers a set of statistical inferences simultaneously or infers a subset of parameters selected based on the observed values. Errors in inference, including confidence intervals that fail to include their corresponding population parameters or hypothesis tests that incorrectly reject the null hypothesis are more likely to occur when one considers the set as a whole.

In empirical studies often dozens or even hundreds of tests are performed for the same data set. When searching for significant test results, one may come up with false discoveries.

Example: $$m$$ different, independent test of significance level $$\alpha>0$$. (Independence means that the test statistics used are mutually independent – this is usually not true in practice). Let’s assume that a common null hypothesis $$H_0$$ holds for each of the $$m$$ tests. Then

$P \begin{pmatrix} \text{Type I error}\\ \text{by at least} \\ \text{one of the m tests} \end{pmatrix} = 1 - (1 - \alpha)^m =: \alpha_m>\alpha$

Therefore, as $$m$$ increases also the probability of a type I error increases:

Number of tests $$m$$ Probability of at least one type I error ($$\alpha_m$$)
1 0.050
3 0.143
5 0.226
10 0.401
100 0.994

Analogous problem: Construction of $$m$$ many $$(1-\alpha)$$ confidence intervals. $P \begin{pmatrix} \text{at least one of the m confidence} \\ \text{intervals does not contain} \\ \text{the true parameter value} \end{pmatrix} = 1 - (1-\alpha)^m>\alpha$

This represents the general problem of multiple comparisons. In practice, it will not be true that all considered test statistics are mutually independent. (This even complicates the problem.) However, we will still have the effect that the probability of at least one falsely significant result increases with the number $$m$$ of tests, but it will not be equal to $$1-(1-\alpha)^m$$.

A statistically rigorous solution of this problem consists in modifying the constructions of tests or confidence intervals in order to arrive at simultaneous tests: $P\begin{pmatrix} \text{Type I error by} \\ \text{at least one of the m tests} \end{pmatrix} \leq \alpha$ or simultaneous confidence intervals: \begin{align*} &P \begin{pmatrix} \text{At least one of the m confidence} \\ \text{intervals does not contain} \\ \text{the true parameter value} \end{pmatrix} \leq \alpha\\[2ex] \Leftrightarrow\quad &P \begin{pmatrix} \text{All confidence intervals} \\ \text{simultaneously contain the} \\ \text{true parameter values} \end{pmatrix} \geq 1 - \alpha \end{align*}

For certain problems (e.g., analysis of variance) there exist specific procedures for constructing simultaneous confidence intervals. However, the only generally applicable procedure seems to be the Bonferroni correction. It is based on Boole’s inequality.

Theorem (Boole): Let $$A_1, A_2, \dots, A_m$$ denote $$m$$ different events. Then $P(A_1 \cup A_2 \cup \dots \cup A_m) \leq \sum_{i=1}^m P(A_i).$ This inequality also implies that: $P(A_1 \cap A_2 \cap \dots \cap A_m) \ge 1 - \sum_{i=1}^m P(\bar A_i),$ where $$\bar A_i$$ denotes the complementary event “not $$A_i$$”.

Example: Bonferroni adjustment for $$m$$ different tests of level $$\alpha^* = \alpha/m$$. $P\begin{pmatrix} \text{Type I error by} \\ \text{at least one of the m tests} \end{pmatrix} \leq \sum_{i=1}^m \alpha^\ast = \alpha$

Analogously: Construction of $$m$$ many $$(1-\alpha^*)$$-confidence intervals with $$\alpha^* =\alpha/m$$: \begin{align*} &P\begin{pmatrix} \text{At least one of the m confidence} \\ \text{intervals does not contain} \\ \text{the true parameter value} \end{pmatrix} \leq \sum_{i=1}^m \alpha^\ast = \alpha\\[2ex] \Leftrightarrow\quad &P\begin{pmatrix} \text{All confidence interval} \\ \text{simultaneously contain the} \\ \text{true parameter values} \end{pmatrix} \geq 1 - \sum_{i=1}^m \alpha^\ast = 1 - \alpha \end{align*}

Example: Regression analysis with $$K=100$$ regressors, where none of the variables has an effect on the dependent variable $$y$$.

library("tidyverse", quietly = TRUE)
K <- 100 # number of regressors
n <- 500 # sample size

set.seed(123)

# Generate regression data, where none of the X-variables
# has an effect on the dependent variable Y:
my_df <- matrix(data = rnorm(n = n*K),
nrow = n, ncol = K,
dimnames = list(paste0("i.",1:n),
paste0("X.",1:K))) %>%
as_tibble() %>%
mutate(Y = rnorm(n)) %>% # Adding a Y-variable that is independent of the X-variables
select(Y, everything())

# OLS regression
OLS_result_df <- lm(Y ~ . , data = my_df) %>%
summary %>%
broom::tidy()

Count_Signif <- OLS_result_df %>%
filter(term != '(Intercept)') %>%
count(p.value < 0.05)
p.value < 0.05 n
FALSE 96
TRUE 4

## 2.6 R-Lab: The Gauss-Test

Let’s reconsider the simplest test statistic you will ever meet: The Gauss-Test (Or “Z-Test”).

Setup: Let $$X_1,\dots,X_n$$ be an iid random sample with $$X_i\sim N(\mu,\sigma^2)$$ and $$\sigma^2<\infty.$$

Idea: Under the above setup, $$\bar{X}_n=n^{-1}\sum_{i=1}^n X_i$$ consistently estimates the (unknown) true mean value $$\mu$$. That is, $$\bar{X}_n\to_p\mu$$.

• Under the null hypothesis (i.e., $$\mu_{0}=\mu$$), the difference $$\bar{X}_n-\mu_{0}$$ should be “small”.
• Under the alternative hypothesis (i.e., $$\mu_{0}\neq\mu$$), the difference $$\bar{X}_n-\mu_{0}$$ should be “large”.

Under the null hypothesis $$H_0$$ we have that $$\mu_{0}=\mu$$. Therefore: $Z=\frac{\sqrt{n}\,(\bar{X}_n-\mu_{0})}{\sigma}=\underbrace{\frac{\sqrt{n}\,(\bar{X}_n-\mu)}{\sigma}}_{\sim N(0,1)}$

Under the alternative $$H_1$$ we have that $$\;\mu_{0}\neq \mu$$. Therefore: \begin{align*} Z&=\frac{\sqrt{n}\,(\bar{X}_n-\mu_{0})}{\sigma}\\ &=\frac{\sqrt{n}\,(\bar{X}_n-\mu_{0}+\mu-\mu)}{\sigma}\\ &=\frac{\sqrt{n}\,(\bar{X}_n-\mu)}{\sigma}+{\color{blue}{ \frac{\sqrt{n}\,(\mu-\mu_{0})}{\sigma} }} \sim N\left({\color{blue}{\frac{\sqrt{n}\,(\mu-\mu_{0})}{\sigma}}},1\right) \end{align*}

The different distributions (under $$H_0$$ and $$H_1$$) of the test statistic $$Z$$ can be investigated in the following dynamic plot:


text:
1.8k views

Consider the following propositional statements:

• $P_1: ((A ∧ B) → C)) ≡ ((A → C) ∧ (B → C))$
• $P_2: ((A ∨ B) → C)) ≡ ((A → C) ∨ (B → C))$

Which one of the following is true?

1. $P_1$ is a tautology, but not $P_2$
2. $P_2$ is a tautology, but not $P_1$
3. $P_1$ and $P_2$ are both tautologies
4. Both $P_1$ and $P_2$ are not tautologies
edited | 1.8k views
+1
sir what does $\equiv$ denotes in this question ??

When i simplify $P1$ ,it comes out to be

$\left ( \left ( A\wedge B \right )\rightarrow C \right )={\left ( AB \right )}'+C$

$A{}'+B{}'+C$ in LHS,

while  In RHS

$\left ( \left ( A\rightarrow C \right ) \wedge \left ( B\rightarrow C \right )\right )$

$\left ( A{}'+C \right ) \wedge \left ( B{}'+C \right )$

=$A{}'B{}'+C$

1.Here LHS$\neq$ RHS ,what is this $\equiv$ here?

2.also i made Kmap of LHS ,i am not getting true(all 1)

while making hte K-Map of LHS ,i am not getting true(all 1) ...
+4
$A\equiv B$ means $\overline{A}$  $\overline{B} + AB$ would be 1.

$\odot$ (EX-NOR)operator is actually 'equivalence operator'
0
A<->B could also be treated as equivalence operator?
0
<-> and ≡ are same I guess.

Also A<->B is same as (A->B)^(B->A).So we need to check whether (A->B)^(B->A) is tautology for P1 and P2.

(D) Both $P_1$ and $P_2$ are not tautologies.

$P_1:$ If $A$ is true and $B$ is false, LHS of $P_1$ is true but RHS becomes false. Hence not tautology.

$P_2:$ Forward side is true. But reverse side is not true. When $A$ is false and $B$ is true and C is false, RHS is true but LHS is false.

LHS of $P_2$ can be simplified as follows:

$((A∨B) → C) \equiv (~(A∨B) ∨ C)$
$\quad \quad \equiv (~A ∧~B) ∨C)$
$\quad \quad \equiv (~A ∨C) ∧ (~B ∨C)$
$\quad \quad \equiv (A→C) ∧ (B→C)$
edited by
0

Can you please covert ((A ∧ B) → C)) to ((A → C) V (B → C))

I can convert into A->(B->C) only

0
they are not equivalent, so how can we convert? You didn't see the last statement- on RHS we get AND in stead of OR, which is a more stricter condition. This makes the forward implication true but reverse implication false.
+1
Sorry..I can't get it.Truth table says they are equivalent.In text book they ask to
d) Transform pq -> r into (p -> r) + (q -> r).
+4
Sorry, I misread the expression. I thought it was P2 in question.
((A ∧ B) → C))  = ~A V ~B V C = ~A V C V ~B V  C = (A → C) V (B → C)
+1
It's ok :) Thank you for the answer
0
if equivalence hold in above, then will it be tautology??
0
sir what does $\equiv$ denotes in this question ??

When i simplify $P1$ ,it comes out to be

$\left ( \left ( A\wedge B \right )\rightarrow C \right )={\left ( AB \right )}'+C$

$A{}'+B{}'+C$ in LHS,

while  In RHS

$\left ( \left ( A\rightarrow C \right ) \wedge \left ( B\rightarrow C \right )\right )$

$\left ( A{}'+C \right ) \wedge \left ( B{}'+C \right )$

=$A{}'B{}'+C$

1.Here LHS$\neq$ RHS ,what is this $\equiv$ here?

2.also i made Kmap of LHS ,i am not getting true(all 1)

while making hte K-Map of LHS ,i am not getting true(all 1) ...
0
I Understand this concept
+1
For proving whether it is a tautology or not...we have to check LHS=RHS or not?? or what else?
0
I can't understand the question,what are they want in the question???

please anyone help me to understand this
0
They want to ask, whether, these expressions are tautologies or not?

$\equiv$ treated as Equivalence
0
Yes.

For equivalence if we prove that LHS=RHS then we can say that proposition is tautology.

Since we consider equicalence as $\odot$ (Ex-nor) we can prove either LHS=RHS=TRUE/FALSE.

$P\leftrightarrow Q$ is same as $P\equiv Q.$ if P=Q=TRUE/FALSE. Then they are equivalent.

Answer is "D" for this question

but what if we exchange in right side of P1 & P2 let's see what happen

P1 : ((A ∧ B) → C)) ≡ ((A → C) ∨ (B → C))
P2 : ((A ∨ B) → C)) ≡ ((A → C) ∧ (B → C))

P1: ~A +~B +C  ≡ ~A +~B+C

P2:  ~A~B + C ≡(~A+C)(~B+C)     // Here we know that C+~A~B = (C+~A)(C+~B)

P1 and P2 are both tautologies here


text:
# Explaining SNARKs Part III: The Knowledge of Coefficient Test and Assumption

<< Part II

In Part II, we saw how Alice can blindly evaluate the hiding :math:E(P(s)) of her polynomial :math:P of degree :math:d, at a point :math:s belonging to Bob. We called this “blind” evaluation, because Alice did not learn :math:s in the process.

However, there was something missing in that protocol – the fact that Alice is able to compute :math:E(P(s)) does not guarantee she will indeed send :math:E(P(s)) to Bob, rather than some completely unrelated value.

Thus, we need a way to “force” Alice to follow the protocol correctly. We will explain in part IV precisely how we achieve this. In this post, we focus on explaining the basic tool needed for that – which we call here the Knowledge of Coefficient (KC) Test.

As before, we denote by :math:g a generator of a group :math:G of order :math:|G|=p where the discrete log is hard. It will be convenient from this post onwards to write our group additively rather than multiplicatively. That is, for :math:\alpha\in\mathbb{F}_p, :math:\alpha\cdot g denotes the result of summing :math:\alpha copies of :math:g.

## The KC Test

For :math:\alpha\in\mathbb{F}_p^* [1], let us call a pair of elements :math:(a,b) in :math:G an :math:\alpha-pair if :math:a,b \neq 0 and :math:b=\alpha\cdot a.

The KC Test proceeds as follows.

1. Bob chooses random :math:\alpha\in\mathbb{F}_p^* and :math:a\in G. He computes :math:b=\alpha\cdot a.
2. He sends to Alice the “challenge” pair :math:(a,b). Note that :math:(a,b) is an :math:\alpha-pair.
3. Alice must now respond with a different pair :math:(a’,b’) that is also an :math:\alpha-pair.
4. Bob accepts Alice’s response only if :math:(a’,b’) is indeed an :math:\alpha-pair. (As he knows :math:\alpha he can check if :math:b’=\alpha\cdot a’.)

Now, let’s think how Alice could successfully respond to the challenge. Let’s assume for a second that she knew :math:\alpha. In that case, she could simply choose any :math:a’ in :math:G, and compute :math:b’=\alpha\cdot a’; and return :math:(a’,b’) as her new :math:\alpha-pair.

However, as the only information about :math:\alpha she has is :math:\alpha\cdot a and :math:G has a hard discrete log problem, we expect that Alice cannot find :math:\alpha.

So how can she successfully respond to the challenge without knowing :math:\alpha?

Here’s the natural way to do it: Alice simply chooses some :math:\gamma\in\mathbb{F}_p^*, and responds with :math:(a’,b’)=(\gamma\cdot a,\gamma\cdot b).

In this case, we have:

:math:b’=\gamma \cdot b = \gamma \alpha \cdot a = \alpha (\gamma\cdot a) =\alpha \cdot a’,

so indeed :math:(a’,b’) is an :math:\alpha-pair as required.

Note that if Alice responds using this strategy, she knows the ratio between :math:a and :math:a’. That is, she knows the coefficient :math:\gamma such that :math:a’=\gamma\cdot a.

The Knowledge of Coefficient Assumption [2] (KCA) states that this is always the case, namely:

KCA: If Alice returns a valid response :math:(a’,b’) to Bob’s challenge :math:(a,b) with non-negligible probability over Bob’s choices of :math:a,\alpha, then she knows :math:\gamma such that :math:a’=\gamma\cdot a.

The KC Test and Assumption will be important tools in Part IV.

## What does “Alice knows” mean exactly

You may wonder how we can phrase the KCA in precise mathematical terms; specifically, how do we formalize the notion that “Alice knows :math:\gamma” in a mathematical definition?

This is done roughly as follows: We say that, in addition to Alice, we have another party which we call Alice’s Extractor. Alice’s Extractor has access to Alice’s inner state.

We then formulate the KCA as saying that: whenever Alice successfully responds with an :math:\alpha-pair :math:(a’,b’), Alice’s Extractor outputs :math:\gamma such that :math:a’=\gamma\cdot a. [3]

 [1] :math:\mathbb{F}_p^* denotes the non-zero elements of :math:\mathbb{F}_p. It is the same as :math:\mathbb{Z}_p^* described in Part I.
 [2] This is typically called the Knowledge of Exponent Assumption in the literature, as traditionally it was used for groups written multiplicatively.
 [3] The fully formal definition needs to give the Extractor “a little slack” and states instead that the probability that Alice responds successfully but the Extractor does not output such :math:\gamma is negligible.

Part IV >>


text:
B. The Feast and the Bus
time limit per test
2 seconds
memory limit per test
512 megabytes
input
standard input
output
standard output

Employees of JebTrains are on their way to celebrate the 256-th day of the year! There are $n$ employees and $k$ teams in JebTrains. Each employee is a member of some (exactly one) team. All teams are numbered from $1$ to $k$. You are given an array of numbers $t_1, t_2, \dots, t_n$ where $t_i$ is the $i$-th employee's team number.

JebTrains is going to rent a single bus to get employees to the feast. The bus will take one or more rides. A bus can pick up an entire team or two entire teams. If three or more teams take a ride together they may start a new project which is considered unacceptable. It's prohibited to split a team, so all members of a team should take the same ride.

It is possible to rent a bus of any capacity $s$. Such a bus can take up to $s$ people on a single ride. The total cost of the rent is equal to $s \cdot r$ burles where $r$ is the number of rides. Note that it's impossible to rent two or more buses.

Help JebTrains to calculate the minimum cost of the rent, required to get all employees to the feast, fulfilling all the conditions above.

Input

The first line contains two integers $n$ and $k$ ($1 \le n \le 5\cdot10^5, 1 \le k \le 8000$) — the number of employees and the number of teams in JebTrains. The second line contains a sequence of integers $t_1, t_2, \dots, t_n$, where $t_i$ ($1 \le t_i \le k$) is the $i$-th employee's team number. Every team contains at least one employee.

Output

Print the minimum cost of the rent.

Examples
Input
6 3
3 1 2 3 2 3

Output
6

Input
10 1
1 1 1 1 1 1 1 1 1 1

Output
10

Input
12 4
1 2 3 1 2 3 4 1 2 1 2 1

Output
12


text:
# Elastic collision of neon atom problem

1. Oct 31, 2008

### mithanon

1. The problem statement, all variables and given/known data
A neon atom (m = 20.0 u) makes a perfectly elastic collision with another atom at rest. After the impact, the neon atom travels away at a 57.9° angle from its original direction and the unknown atom travels away at a -45.4° angle. What is the mass (in u) of the unknown atom? [Hint: You can use the law of sines.]

variables
m = neon = 20 u
M = mass of unknown
$$\theta$$ = 57.9°
$$\phi$$ = 45.4°
v = initial velocity of neon
v' = final velocity of neon
V = final velocity of unknown

2. Relevant equations
I've set up conservation of momentum equations for both directions, and the energy equation, but I'm stuck with more unknowns than equations

3. The attempt at a solution
mv2=mv'2+MV2
mv = mv'cos$$\theta$$ + MVcos$$\phi$$
0 = mv'sin$$\theta$$ - MVsin$$\phi$$

using a momentum vector sum diagram and law of sines? I got
mv/sin(180-$$\theta$$-$$\phi$$)=mv'/sin($$\phi$$)

not sure how that helps but I guess I can then express v' in terms of v, eliminating one variable

no idea where to go from here

2. Nov 1, 2008

### hage567

Everything looks OK so far I think. You can also use the law of sines to get the other "side" of the triangle. The problem is now just one of algebra. I would pick one equation (maybe the energy one), and substitute in for the variables you don't want and then isolate M.


text:
# Triangulation with maximum greatest area

Given a point set $P$ and a triangulation $T$ of $P$ with $d$ triangles, let's define

$$\alpha(T) = (\alpha_1, \alpha_2, \ldots, \alpha_{3d})$$

which denotes the series of interior angles of $T$, sorted in ascending order. We can compare two such series (of the same length) easily in lexicographical order. We know that $T$ is a Delaunay-Triangulation iff $\alpha(T) \geq \alpha(T')$ for all possible triangulations $T'$ of $P$. This is also called the "maximum smallest angle" characteristic of the Delaunay Triangulation.

My question is now, is there a "maximum greatest area" triangulation? (Yes, not "maximum smallest area"!) More formally, define

$$A(T) = (A_1, A_2, \ldots, A_d)$$

which denotes the series of triangle areas of $T$, sorted in descending order, and I hope to find a triangulation $T$, so that $A(T) \geq A(T')$ for all possible triangulations $T'$ of $P$.

Does anyone know such an algorithm, or is it even trivially computable (i.e. with a greedy algorithm)?

• I think as you suggest it's trivially computable by a greedy algorithm when there are no ties ub the area between two different triangles. When there are ties, it might be more difficult. – David Eppstein Mar 25 '12 at 16:43


text:
How to Break Coupled ODEs down to first order for Runge-Kutta

My question might seem a bit simple. I am trying to solve a system of ODEs using Runge-Kutta method. I am having difficulty breaking down the equations into a system of first order ones required before applying RK45 to it because there is a third order differential of the same parameter "F" in both equations which I am not used to.

Equation 1 $$\left( {\frac {{\rm d}^{2}}{{\rm d}{\eta}^{2}}}F \left( \eta \right) \right) ^{n-1}{\frac {{\rm d}^{3}}{{\rm d}{\eta}^{3}}}F \left( \eta \right) +{\frac {2\,n+2+ \left( 2\,n-1 \right) kF \left( \eta \right) {\frac {{\rm d}^{2}}{{\rm d}{\eta}^{2}}}F \left( \eta \right) }{n+4}}-{\frac { \left( n+2+ \left( n+1 \right) k \right) \left( {\frac {\rm d}{{\rm d}\eta}}F \left( \eta \right) \right) ^{2 }}{n+4}}+G \left( \eta \right) =0$$

Equation 2

$${\frac { \left( {\frac {{\rm d}^{2}}{{\rm d}{\eta}^{2}}}F \left( \eta \right) \right) ^{n-1} \left( n-1 \right) \left( {\frac {{\rm d}^{3 }}{{\rm d}{\eta}^{3}}}F \left( \eta \right) \right) {\frac {\rm d}{ {\rm d}\eta}}G \left( \eta \right) }{\Pr\,{\frac {{\rm d}^{2}}{{\rm d} {\eta}^{2}}}F \left( \eta \right) }}+ \left( {\frac {{\rm d}^{2}}{ {\rm d}{\eta}^{2}}}F \left( \eta \right) \right) ^{n-1}{\frac { {\rm d}^{2}}{{\rm d}{\eta}^{2}}}G \left( \eta \right) +{\frac {2\,n+2+ \left( 2\,n-1 \right) kF \left( \eta \right) {\frac {\rm d}{{\rm d} \eta}}G \left( \eta \right) }{n+4}}-{\frac {n+2\, \left( n+1 \right) k \left( {\frac {\rm d}{{\rm d}\eta}}F \left( \eta \right) \right) G \left( \eta \right) }{n+4}}=0$$

Boundary Conditions

$$G \left( 0 \right) =1, F \left( 0 \right)=F' \left( 0 \right) =0,G \left( \infty \right) =0,F' \left( \infty \right) =0$$

Any help would be greatly appreciated.

• Out of curiosity, is there any physical meaning behind these equations? They look pretty arbitrary to me. I tried several times to change variables and make it more simpler but it fails with no luck... There are some patterns in these two equations that encouraged me to try make it simpler or write it down as a system of first-order ODEs but at least I can't get it... Are you sure the constants in these two equations are written correctly? For example you have: $\frac{n+2+(n+1)k}{n+4}$ and $\frac{n+2(n+1)k}{n+4}$, which looks pretty similar but not exactly the same... – Alone Programmer Oct 3 '19 at 21:24
• If you are using RK45, you are not going to be able to enforce the boundary conditions at $\infty$. You are going to have to get lucky with your choice of initial data – whpowell96 Oct 3 '19 at 23:43

The usual trick is to add more variables that represent the successive derivatives, as in the equations of motion pf physics written as a set of first order ODE of "2" variables: \begin{align}\dot{x}&=v\\ \dot{v}&=a(x,t)\end{align} instead of a second order ODE of 1 variable $$\ddot{x}=a(x,t) \ .$$ The dot represents the time derivative.

So, applying to your case, you could have Equation 1 as: \begin{align} x & = \frac{d F}{d \eta}\\ y & = \frac{d x}{d \eta} \\ 0 & = y^{n-1}\frac {d y}{d \eta} + \frac{2n+2 + (2n-1) k F y}{n+4} - \frac { (n+2 + (n+1) k) x^2}{n+4} + G(\eta) \end{align} Note that $$y = \frac{d x}{d \eta} = \frac{d^2 F(\eta)}{d \eta^2}$$ and that the last equation can easily be written as $$\frac{d y}{d \eta} = \ldots$$

Equation 2 can have similar substitutions together with the previous ones: \begin{align} u & = \frac{d G}{d \eta}\\ v & = \frac{d u}{d \eta} \\ 0 & = \frac{y^{n-1} (n-1) \frac{d y}{d \eta} u }{\Pr\,y} + y^{n-1} v + \frac {2n+2 + (2n-1) k F u}{n+4} - \frac {n+2 (n+1) k x G}{n+4} \end{align} And the last equation can also be written as $$\frac{d y}{d \eta} = \ldots$$

This means that you will have 4 virtual new variables: $$x$$, $$y$$, $$u$$ and $$v$$.

Edit:

Looking at the number of boundary conditions, since you have 5, you might want to have 5 variables. Depending on your problem, you may write Equation 2 as: \begin{align} u & = \frac{d G}{d \eta}\\ 0 & = \frac{y^{n-1} (n-1) \frac{d y}{d \eta} u }{\Pr\,y} + y^{n-1} \frac{d u}{d \eta} + \frac {2n+2 + (2n-1) k F u}{n+4} - \frac {n+2 (n+1) k x G}{n+4} \end{align} The difference is in the second parcel, where we have $$\frac{d u}{d \eta}$$ and we no longer need to have $$v$$.

And then, plug in the expression you got for $$\frac{d y}{d \eta}$$ in the last equation and solve it for $$\frac{d u}{d \eta} = \ldots$$

Hence, you'll get a system with 5 (first order) ODE and 5 boundary conditions.