Advanced Search

Alec Mihailovs

Dr. Aleksandrs Mihailovs

4495 Reputation

21 Badges

20 years, 335 days
Mihailovs, Inc.
Owner, President, and CEO
Tyngsboro, Massachusetts, United States

Social Networks and Content at Maplesoft.com

Maple Application Center
Contact Dr. Aleksandrs Mihailovs

I received my Ph.D. from the University of Pennsylvania in 1998 and I have been teaching since then at SUNY Oneonta for 1 year, at Shepherd University for 5 years, at Tennessee Tech for 2 years, at Lane College for 1 year, and this year I taught at the University of Massachusetts Lowell. My research interests include Representation Theory and Combinatorics.

MaplePrimes Activity


These are replies submitted by Alec Mihailovs

Special cases...

Alec Mihailovs 4495
September 15 2005
Yes, in general, Cholesky decomposition gives a solution. In this example, with c[i,j] = rho for i ≠ j, it is rather cumbersome. I used a different approach. If there exist such numbers c[i] that c[i,j] =c[i]*c[j] for i ≠ j, then the samples B[i] with correlation matrix c can be constructed as 
               B[i] = c[i]* W[0] + sqrt(1-c[i]^2)* W[i]
In this example, c[i] = sqrt(rho).

Works OK...

Alec Mihailovs 4495
September 15 2005
I just tried a worksheet that I haven't opened earlier, so it is not in the cache, and it worked OK.

Right...

Alec Mihailovs 4495
September 15 2005
0 0
Yes, it is the right way to do that. It can't be compared with evalhf/zip in previous versions, because neither Statistics nor evalhf/zip existed in previous versions. Here are time comparisons for Axel Vogt's example, 
B:=[stats[random,normald[0,1]](10^5)]:
W:=[stats[random,normald[0,1]](10^5)]:

time(zip((x,y) -> evalf(rho*x+sqrt(1-rho^2)*y) ,B,W)); 

                                22.983

time(zip((x,y) -> evalhf(rho*x+sqrt(1-rho^2)*y) ,B,W)); 

                                2.718

time(rho*B+sqrt(1-rho^2)*W);

                                0.014

Right...

Alec Mihailovs 4495
September 15 2005
Yes, it is the right way to do that. It can't be compared with evalhf/zip in previous versions, because neither Statistics nor evalhf/zip existed in previous versions. Here are time comparisons for Axel Vogt's example, 
B:=[stats[random,normald[0,1]](10^5)]:
W:=[stats[random,normald[0,1]](10^5)]:

time(zip((x,y) -> evalf(rho*x+sqrt(1-rho^2)*y) ,B,W)); 

                                22.983

time(zip((x,y) -> evalhf(rho*x+sqrt(1-rho^2)*y) ,B,W)); 

                                2.718

time(rho*B+sqrt(1-rho^2)*W);

                                0.014

100 datasets with given correlation...

Alec Mihailovs 4495
September 15 2005
0 0
If ρ ≥ 0, one can obtain 100 datasets B1, ..., B100 with mutual correlation ρ as follows, 
with(Statistics): 
X := RandomVariable(Normal(0, 1)): 
rho:=0.5:
A := Sample(X, 10^5):
for n to 100 do B[n]:=sqrt(rho)*A+sqrt(1-rho)*Sample(X,10^5) od:
Test it, 
s:=seq(seq(Correlation(B[i],B[j]),i=j+1..100),j=1..99):
min(s),max(s);
                      0.4773409752, 0.5261740852,
With 100 values, the correlations may be not that close to ρ. For example, 
A := Sample(X, 100):
for n to 100 do B[n]:=sqrt(rho)*A+sqrt(1-rho)*Sample(X,100) od:
s:=seq(seq(Correlation(B[i],B[j]),i=j+1..100),j=1..99):
min(s),max(s);
                      0.2237998717, 0.6950226961

100 datasets with given correlation...

Alec Mihailovs 4495
September 15 2005
If ρ ≥ 0, one can obtain 100 datasets B1, ..., B100 with mutual correlation ρ as follows, 
with(Statistics): 
X := RandomVariable(Normal(0, 1)): 
rho:=0.5:
A := Sample(X, 10^5):
for n to 100 do B[n]:=sqrt(rho)*A+sqrt(1-rho)*Sample(X,10^5) od:
Test it, 
s:=seq(seq(Correlation(B[i],B[j]),i=j+1..100),j=1..99):
min(s),max(s);
                      0.4773409752, 0.5261740852,
With 100 values, the correlations may be not that close to ρ. For example, 
A := Sample(X, 100):
for n to 100 do B[n]:=sqrt(rho)*A+sqrt(1-rho)*Sample(X,100) od:
s:=seq(seq(Correlation(B[i],B[j]),i=j+1..100),j=1..99):
min(s),max(s);
                      0.2237998717, 0.6950226961

evalhf(zip(......

Alec Mihailovs 4495
September 14 2005
zip may be slow, but evalhf(zip(... is quite fast.

Faster...

Alec Mihailovs 4495
September 14 2005
0 0
It is faster with new Statistics package, 
with(Statistics): 
X := RandomVariable(Normal(0, 1)): 
Y := RandomVariable(Normal(0, 1)): 
rho:=0.5:
A := Sample(X, 10^5):
B := Sample(Y, 10^5):
C := evalhf(zip((x,y)->rho*x+sqrt(1-rho^2)*y,A,B)):

Faster...

Alec Mihailovs 4495
September 14 2005
It is faster with new Statistics package, 
with(Statistics): 
X := RandomVariable(Normal(0, 1)): 
Y := RandomVariable(Normal(0, 1)): 
rho:=0.5:
A := Sample(X, 10^5):
B := Sample(Y, 10^5):
C := evalhf(zip((x,y)->rho*x+sqrt(1-rho^2)*y,A,B)):

Goto...

Alec Mihailovs 4495
September 13 2005
In addition to undocumented library procedures, there are quite a few undocumented built-in commands, such as goto. For example, 
a:=proc()
goto(1);
print("not interesting");
1: print("interesting") end:

> a();
                            "interesting"
Joe Riel is the best expert in the undocumented Maple features, so his comment on this topic would be the most appropriate.

A solution...

Alec Mihailovs 4495
September 12 2005
0 0
The numerical solution can be found in a similar way, just with interchanging of x0 and 0. 
gfx0:=unapply(solve(ug - f0 = c2*sqrt( g(f0) - g(f(x0)) ),g(f(x0))),f0):
fx0:=f0->rhs(dsolve({diff(f(x),x) = c1*sqrt( g(f(x)) - gfx0(f0)), f(0)=f0},
numeric,method=dverk78)(x0)[2]):
Now it is time to enter the parameters, 
g := f->exp(f-xn) +exp(-f) +(1-exp(-xn))*f:
c1,c2,xn,ug,x0:=1,-2,5,1,-1:
plot({g@fx0,gfx0},1.8..2.2);
f0:=fsolve(g('fx0'(f0))-gfx0(f0),f0=1.85);

                          f0 := 1.849058144

sol:=dsolve({diff(f(x),x) = c1*sqrt( g(f(x)) - gfx0(f0)), f(0)=f0},
numeric,method=dverk78):
plot(x->rhs(sol(x)[2]),-2..4,0..9);

A solution...

Alec Mihailovs 4495
September 12 2005
The numerical solution can be found in a similar way, just with interchanging of x0 and 0. 
gfx0:=unapply(solve(ug - f0 = c2*sqrt( g(f0) - g(f(x0)) ),g(f(x0))),f0):
fx0:=f0->rhs(dsolve({diff(f(x),x) = c1*sqrt( g(f(x)) - gfx0(f0)), f(0)=f0},
numeric,method=dverk78)(x0)[2]):
Now it is time to enter the parameters, 
g := f->exp(f-xn) +exp(-f) +(1-exp(-xn))*f:
c1,c2,xn,ug,x0:=1,-2,5,1,-1:
plot({g@fx0,gfx0},1.8..2.2);
f0:=fsolve(g('fx0'(f0))-gfx0(f0),f0=1.85);

                          f0 := 1.849058144

sol:=dsolve({diff(f(x),x) = c1*sqrt( g(f(x)) - gfx0(f0)), f(0)=f0},
numeric,method=dverk78):
plot(x->rhs(sol(x)[2]),-2..4,0..9);

Solutions...

Alec Mihailovs 4495
September 12 2005
0 0
Well, IVP {diff(f(x),x) = c1*sqrt( g(f(x)) - g(f0) ), f(x0)=f0} has constant solution f(x)=f0. Do you think it has other solutions?

Solutions...

Alec Mihailovs 4495
September 12 2005
Well, IVP {diff(f(x),x) = c1*sqrt( g(f(x)) - g(f0) ), f(x0)=f0} has constant solution f(x)=f0. Do you think it has other solutions?

Constant solution...

Alec Mihailovs 4495
September 12 2005
0 0
The system of equations DEq and Eq has constant solution f(x)=ug. Here is the description of posting worksheets.
First 169 170 171 172 173 174 175 Last Page 171 of 180