[R] Integrate and mapply

[Joshua Wiley]

Hi,

I still do not really understand, but I have a couple ideas:

On Sun, Nov 7, 2010 at 11:09 AM, Vaiva P <asvaiva at gmail.com> wrote:
> Thank you guys for  a quick reaction. I decided it to write a common
> letter to you both at once.
> I tried to calculate the value of a a function without integration. I
> gave the value for u and I got nice 15 dimensions vector. That is the
> thing, I need to get with an integral. However, as soon as put
> integral in, I get various error messages.
> One thing that I tried, I tried to put function undint<-function(u)
> inside function prob<-function(i), i.e.
>
> prob<-function(i) {
> i1<-comb[i,1]
> ------------
>  undint<-function(u)
>
> Then the last line of function prob(i) would be
> integrate(unind, 0, Inf)$value }
>
> Afterwards I tried to apply sapply()  function just with the vector of
> i values for function prob(i), but it did not work.
> Error message:
> There were 15 warnings (use warnings() to see them)
> Warning messages:
> 1: In matrix(u, length(val1), 1) :
>  data length [15] is not a sub-multiple or multiple of the number of rows [2]

This sounds like you tried to integrate using prob(), but that is not
a vectorized function.  Consider:

myfun <- function(x) {print(x); return(x+2)}
integrate(f = myfun, lower = 0, upper = 4)

integrate() passes a vector of values at once to myfun() (which you
can see because myfun prints them to the screen).  In your case, it
looks like 15 values were passed to the "u" argument of your function.
 So, if you want to get your prob() function working with integrate,
you first need to be able to have it work when u is a vector.  For
example:

u <- 1
prob(i = 1) # works
u <- seq(0, 4, .4)
prob(i = 1) # fails right now

It should not be terribly difficult to do.  A couple notes that might
help get you started: Vectorize() will not help, because you need
prob() to be able to handle a vector "u", not "i" (its only argument).
 Making "i" a vector will not solve your problem for the same reason.

HTH,

Josh

>
>
> As soon as I delete integrate() and give a numerical value for u, I
> get  nice numerical vector, but I cannot get it if there is an
> integral.
> This integral is  the part of Discrete choice model estimation code.
> In fact I need to define my system of equations that I use with
> nlseqslv() later. I defined 'nu' as some number in the example, not to
> put all clumsy code. In fact it is the unknown that needs to be solved
> for (conceptually similar to 'delta' in BLP model)
>  I tried to make the problem description as simple as possible in my
> letter, because it shows the problem, that I detected after trying all
> parts of the code.
> If some more details are needed I am ready to give them.
>
> Thanks,
> Vaiva
>
> 2010/11/7 Uwe Ligges <ligges at statistik.tu-dortmund.de>:
>> So undint(u) is a 15 dimensional vector. What do the different dimensions
>> mean? How would you define the integral of a 15 dimensional vector? It would
>> help if you could provide some background on what your code is supposed to
>> do.
>>
>> Uwe Ligges
>>
>>
>>
>>
>>
>> On 07.11.2010 17:01, Vaiva P wrote:
>>>
>>> Hi,
>>>
>>> I need some help on integrating a function that is a vector.
>>> I have a function - vector which each element is different. And,
>>> naturally, function integrate() does not work
>>> I checked the article of U. Ligges and J. Fox (2008) about code
>>> optimization "How Can I Avoid This Loop or Make It Faster?" on
>>> http://promberger.info/files/rnews-vectorvsloops2008.pdf.
>>> Their advice did not help me, because they advised a solution, when
>>> function was the same, but limits of the integral were  vectors. I
>>> have got everything other way around. I could not get any result.
>>> I would be very grateful for an advice on this matter.
>>> The code is below. The function needed to be integrated (undint with
>>> respect to u from 0 to Inf).
>>>
>>> Sincerely,
>>> V.
>>>
>>> Code
>>>
>>> rm=list(ls())
>>>
>>> #index matrix
>>> t1<-c(1:15)
>>> t2<-c(1:5)%x%matrix(1,3,1)
>>> t3<-matrix(1,3,1)%x%c(1:5)
>>> t4<-c(1:3)%x%matrix(1,5,1)
>>> comb<-cbind(t1,t2,t3,t4)
>>>
>>> nu<-rchisq(15,4)
>>>
>>> gam<-matrix(rchisq(75,df=5),15,5)
>>>
>>> undint<-function(u){
>>> prob<-function(i) {
>>> i1<-comb[i,1]
>>> i2<-comb[i,2]
>>> i3<-comb[i,3]
>>> i4<-comb[i,4]
>>>
>>> val1<-gam[((i3-1)*2+1):(i3*2),i4]
>>> vc1<-matrix(12,length(val1),1)-matrix(u,length(val1),1)-val1
>>> distr1<-prod(pchisq(vc1,df=2))
>>>
>>>
>>> p10<-distr1*1/(exp(-nu[i1])*(1+gam[i2,i4]))*(1-plogis((u+log(gam[i3,i4])-log(1+gam[i2,i4])),0,1))
>>> p10
>>> }
>>> val<-c(1:15)
>>> q<-sapply(val,prob)
>>> }
>>>
>>> ______________________________________________
>>> R-help at r-project.org mailing list
>>> https://stat.ethz.ch/mailman/listinfo/r-help
>>> PLEASE do read the posting guide
>>> http://www.R-project.org/posting-guide.html
>>> and provide commented, minimal, self-contained, reproducible code.
>>
>



-- 
Joshua Wiley
Ph.D. Student, Health Psychology
University of California, Los Angeles
http://www.joshuawiley.com/