Path: news.mathworks.com!not-for-mail
From: <HIDDEN>
Newsgroups: comp.soft-sys.matlab
Subject: Re: How to avoid for loops in 3 vector function.
Date: Wed, 4 May 2011 11:10:05 +0000 (UTC)
Organization: California Institute of Technology
Lines: 96
Message-ID: <iprc6c$67r$1@fred.mathworks.com>
References: <ipmluu$l21$1@fred.mathworks.com> <ipph4t$8g9$1@fred.mathworks.com> <ippm7h$8au$1@fred.mathworks.com>
Reply-To: <HIDDEN>
NNTP-Posting-Host: www-01-blr.mathworks.com
Content-Type: text/plain; charset=UTF-8; format=flowed
Content-Transfer-Encoding: 8bit
X-Trace: fred.mathworks.com 1304507405 6395 172.30.248.46 (4 May 2011 11:10:05 GMT)
X-Complaints-To: news@mathworks.com
NNTP-Posting-Date: Wed, 4 May 2011 11:10:05 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 2843808
Xref: news.mathworks.com comp.soft-sys.matlab:725158

"Roger Stafford" wrote in message <ippm7h$8au$1@fred.mathworks.com>...
> "Roger Stafford" wrote in message <ipph4t$8g9$1@fred.mathworks.com>...
> > "Lan" wrote in message <ipmluu$l21$1@fred.mathworks.com>...
> > > I have a problem with a function Sn(M0, V0, v, b). M0, V0 are 2000x1 given vectors. I want to maximize Sn over b, for each different v belonging to V0.
> > > 
> > > V0, M0, Y are n x 1 vectors (n=2000)
> > > 
> > > for i=1:length(V0)
> > >     v=V0(i);
> > >     b = (-3:0.01:3);
> > >     l = length(b);
> > >     vec0 = 2*Y - ones(n,1);% n x 1 vector
> > >     mV0=repmat(V0,1,l); % n x n matrix
> > >     t=repmat(M0, 1, l)+v*ones(n, l)+repmat(b,n,1); % n x l matrix
> > >     vec1=(t./h+0.5*ones(n,l)).*(abs(t./h) <= 0.5) +((t./h) > 0.5); % n x l matrix
> > >     vec2=(abs((mV0 - v*ones(n,l))/h)<=0.5); % n x l matrix
> > >     sn0=vec0'*(vec1.*vec2);
> > >     sn=1/(n*h)*sn0; % 1 x l array
> > >     p = find(sn>=max(sn));
> > >     betaH(i)=b(p(1));
> > > end
> > > 
> > > This runs very slowly. How do I vectorize the loop? Or is there any other way to make the function run faster?
> > - - - - - - - - - - -
> >   You have defined vector b as the 601 discrete elements of the vector -3:0.01:3, which makes me think that in the ideal model you are working on, you would very likely prefer to consider the entire infinitude of all points within the interval -3<=b<=3 but have settled for the discrete case only because of practical computational considerations.  If that guess is correct, it should be possible to revert back to this ideal continuum model and at the same time greatly increase the efficiency of your computation, by working with interval bounds in b rather than predefined discrete values.
> > .......
> - - - - - - - -
>   I withdraw that suggestion I made earlier in this thread.  In my haste I misinterpreted the vec2 vector as being entirely logical with only 0 and 1 values.  Also I didn't interpret sn0 correctly.  I am wondering if your 'h' is somehow related to the intervals between values in vector 'b'.
> 
>   I add my voice to Florin's and urge you to describe your problem to us in more descriptive terminology if we are to effectively help you.  It is difficult to get a good grasp on the problem from the abstract quantities you have defined.
> 
> Roger Stafford


Thanks for all advice. I will try to describe the whole problem as clearly as possible. 

This is the problem of simulating the binary response model. My binary response model is Y(j) = ((x0(j)+b0(j))>0). Remark: this is a scalar example , for the future, I would like to model Y(j)= (x(j)*b(j)>0) ,0 <j<=n, where x(j), and b(j) are vectors of dimension d. and * is a scalar product. But for now let's just focus on the concrete scalar example.

I am first generating multivariate random variables b0, z0, and x0. 
x0 is endogenous, z0 is an instrument. The goal is to estimate b from the sample, by maximizing Sn function, defined later on, locally. 

mu = [0 0 0];
SIGMA = [1 0.7 0; 0.7 1 0.6; 0 0.6 1];
n=2000;
M = mvnrnd(mu,SIGMA,n);
b0 = M(:,1);
x0 = M(:,2);
z0 = M(:,3);

We have a relation: x0 = M0(z) + V0(x, z). z is not correlated with b.

I compute M0

for n=1:length(z0)
    X=-100:2/100:100;
    Y=X.*normpdf(X,(0.6*z0(n)),0.64);
    M0(n) = trapz(X,Y);
end

V0 = x0-M0;

Y = (x0+b0 > 0);
betaH = zeros(n,1);

b = (-3:0.1:3);
l = length(b);
vec0 = 2*Y - ones(n,1);% n x 1 array
mV0=repmat(V0,1,l); % n x n matrix
t=zeros(n,l);
vec1=zeros(n,l);
vec2=zeros(n,l);
t_temp=repmat(M0, 1, l)+repmat(b,n,1);
my_ones=ones(n,l);


for i=1:length(V0)
    v=V0(i);
    t=t_temp+v*my_ones;
    % t=repmat(M0, 1, l)+v*ones(n, l)+repmat(b,n,1); % n x l matrix
    vec1=(t./h+0.5*ones(n,l)).*(abs(t./h) <= 0.5) +((t./h) > 0.5); % n x l matrix
    vec2=(abs((mV0 - v*ones(n,l))/h)<=0.5); % n x l matrix
    sn0=vec0'*(vec1.*vec2);
    Sn=1/(n*h)*sn0; % 1 x l array
    p = find(Sn>=max(sn));
    betaH(i)=b(p(1));
end


vec2 represents a matrix of integrated kernel function evaluated at every b.
vec1 is a kernel function.

Maximizing Sn(v,b) over b for a given v, for every v, will yield a relation betaH(v). averaging the calculated betaH will give an estimate of 'true' b. For Now, I would like to be able to plot betaH vs V0.

I hope the above made sense somehow. Let me know if something is confusing.

Lan