Create a function with multiple parameters in R -

i want compute following functions :

here, g(x) density function of distribution. want compute function several distributions. in addition, use library fitdistrplus.

to create g, use function do.call way :

g<-function(x) {do.call(paste("d",i,sep=""),c(list(x=x),fti$estimate))} 

fti$estimate contains parameters of distribution i.

g(x) cumulative distribution computed way :

g<-function(x) {do.call(paste("p",i,sep=""),c(list(q=x),fti$estimate))} 

i compute f(x) way :

f<function(n,x) {n*g(x)*(1-g(x))^(n-1) 

at last, compute h(x) way :

h<- function(n) {integrate(function(x) {x*f(n,x)},0,inf)} 

however, can't plot these functions, following errors :

1: in n*g(x): longer object length not multiple of shorter object length 2: in (1-g(x))^(n-1): longer object length not multiple of shorter object length 3: in x*f(n,x) : longer object length not multiple of shorter object length 

beyond, if juste want plot f(n,x), error :

error in list(x=x) :'x' missing 

the minimal snipset have following

#i can "exp" "lnorm" "norm" etc... for( in functionsname) {     png(paste(filebase,"_",i,"_","graphics.png",sep=""))     plot.new()      fti<-fitdist(data, i)     plotdist(data,i, para=as.list(fti[[1]]))     #fti datatable or datafram     #fti$estimate looks :     # meanlog    sdlog      #8.475449 1.204958       #g     pdf<-function(x) {do.call(paste("d",i,sep=""), c(list(x=x),fti$estimate))}   #g     cdf<-function(x) do.call(paste("p",i,sep=""), c(list(q=x),fti$estimate))   #f     minlaw<- function(n,x) {n*pdf(x)*(1-cdf(x))^(n-1)}     #h     minexpectedvalue<-function(n) {integrate(function(x) {x*minlaw(n,x)},0,inf)}      #these 2 following lines give error     plot(minexpectedvalue)     plot(minlaw)      dev.off() } 

i had reverse engineering figure out d1, q1 etc calls, think how it. perhaps original problem lies in function call f(n=2:3, x=1:9); in such call n should single value, not vector of values.

even if length of x multiple of n length, output not wanted. if try give n vector form, might end in recycled (false) output:

> print(data.frame(n=2:3, x=1:6)) - n x 1 2 1 2 3 2 3 2 3 4 3 4 5 2 5 6 3 6 

where x evaluated n=2 @ point x=1, n=3 @ point x=2 etc. would've wanted in lines of:

> print(expand.grid(x=1:5, n=2:3)) -  x n 1  1 2 2  2 2 3  3 2 4  4 2 5  5 2 6  1 3 7  2 3 8  3 3 9  4 3 10 5 3 

you calling separately each n value:

lapply(2:3, fun=function(n) (f(n, x=1:5))) #[[1]] #[1] 0.0004981910 0.0006066275 0.0007328627 0.0008786344 0.0010456478 # #[[2]] #[1] 0.0007464956 0.0009087272 0.0010974595 0.0013152213 0.0015644676 

did use same fti distribution fits, though should've been different? or i in fti refer index i , list of fits in form of ft[[i]]?

below wrapper function, called separately each n-value (and distribution i):

wrapper <- function(i, x, n, fti){     # provided op     g<-function(x) {do.call(paste("d",i,sep=""),c(list(x=x),fti$estimate))}      g<-function(x) {do.call(paste("p",i,sep=""),c(list(q=x),fti$estimate))}     # in fti refer fit of i:th distribution, i.e. should list i:th location in ft i:th distribution estimates?      f<-function(n,x) {n*g(x)*(1-g(x))^(n-1)}     # missing '-' , '}'      h<- function(n) {integrate(function(x) {x*f(n,x)},0,inf)}      list(gres = g(x), gres = g(x), fres = f(n,x), hres = h(n)) }  # example data require("fitdistrplus") data(groundbeef) serving <- groundbeef$serving  # gumbel distribution d1 <- function(x, a, b) 1/b*exp((a-x)/b)*exp(-exp((a-x)/b)) p1 <- function(q, a, b) exp(-exp((a-q)/b)) q1 <- function(p, a, b) a-b*log(-log(p))  fti1 <- fitdist(serving, "1", start=list(a=10, b=10)) #> fti1$estimate #              b  #56.95893 29.07871  # normal distribution  # dnorm, pnorm , qnorm available in default environment d2 <- dnorm p2 <- pnorm q2 <- qnorm  fti2 <- fitdist(serving, "2", start=list(mean=0, sd=1)) #> fti2$estimate #    mean       sd  #73.67743 35.92581  # sequence of x-values xs <- seq(-100, 100, by=1)  print((resultdist1n2 <- wrapper(i=1, x=xs, n=2, fti=fti1))$hres) print((resultdist1n3 <- wrapper(i=1, x=xs, n=3, fti=fti1))$hres) print((resultdist2n2 <- wrapper(i=2, x=xs, n=2, fti=fti2))$hres) print((resultdist2n3 <- wrapper(i=2, x=xs, n=3, fti=fti2))$hres)  plot(xs, resultdist1n2$fres, col=1, type="l", ylim=c(0,0.025), xlab="x", ylab="f(n, x)") points(xs, resultdist1n3$fres, col=2, type="l") points(xs, resultdist2n2$fres, col=3, type="l") points(xs, resultdist2n3$fres, col=4, type="l") legend("topleft", legend=c("gamma (i=1) n=2", "gamma (i=1) n=3", "normal (i=2) n=2", "normal (i=2) n=3"), col=1:4, lty=1) 

and results of desired h found in resultdist1n2$hres etc:

h(n=2) distribution i=1: 53.59385 absolute error < 0.00022 h(n=3) distribution i=1: 45.23146 absolute error < 4.5e-05 h(n=2) distribution i=2: 53.93748 absolute error < 1.1e-05 h(n=3) distribution i=2: 44.06331 absolute error < 2e-05 

edit: here's how 1 uses lapply function call each of vector of n values 0<=n<=256:

ns <- 0:256 res1 <- lapply(ns, fun=function(nseq) wrapper(i=1, x=xs, n=nseq, fti=fti1)) par(mfrow=c(1,2)) plot.new() plot.window(xlim=c(-100,100), ylim=c(0, 0.05)) box(); axis(1); axis(2); title(xlab="x", ylab="f(n,x)", main="f(n,x) gamma (i=1), n=0:256") for(i in 1:length(ns)) points(xs, res1[[i]]$fres, col=rainbow(257)[i], type="l") # perform other distributions calling i=2, fti=fti2 # h function of n dist i=1 plot(ns, unlist(lapply(res1, fun=function(x) x$hres$value)), col=rainbow(257), xlab="n", ylab="h(n)", main="h(n) gamma (i=1), n=0:256") 

i plot each distribution i separately this.