> | with(linalg); |
> | M:=matrix(3,4,[60,50,53,28,110,120,127,59,30,30,30,15]); |
> | gausselim(M); |
> | gaussjord(M); |
> | P:=Vector[row]([5/10,1/10,4/10]); |
> | m:=[(6/25)/P[1],(3/50)/P[2],(1/5)/P[3]]; |
> | R1:=Vector[row]([M[1,1]/M[1,4],M[1,2]/M[1,4],M[1,3]/M[1,4]]); |
> | R2:=Vector[row]([M[2,1]/M[2,4],M[2,2]/M[2,4],M[2,3]/M[2,4]]); |
> | R3:=Vector[row]([M[3,1]/M[3,4],M[3,2]/M[3,4],M[3,3]/M[3,4]]); |
> | ER1:=dotprod(R1,P);#E(R1) |
> | evalf(%); |
> | ER2:=dotprod(R2,P);#E(R2) |
> | evalf(%); |
> | dotprod(R3,P); |
> | s1:=sqrt(P[1]*(R1[1]-ER1)^2+(P[2])*(R1[2]-ER1)^2+(P[3])*(R1[3]-ER1)^2);#stdev(R1) |
> | evalf(%); |
> | s2:=sqrt(P[1]*(R2[1]-ER2)^2+(P[2])*(R2[2]-ER2)^2+(P[3])*(R2[3]-ER2)^2);#stdev(R2) |
> | evalf(%); |
> | R1R2:=[R1[1]*R2[1],R1[2]*R2[2],R1[3]*R2[3]]; |
> | ER1R2:=dotprod(R1R2,P);#E(R1R2) |
> | (ER1R2-ER1*ER2)/(s1*s2);#rho(R1,R2) |
> | evalf(%); |
> | mR1:=[m[1]*R1[1],m[2]*R1[2],m[3]*R1[3]]; |
> | dotprod(mR1,P);#check that E(mR1)=1 |
> | mR2:=[m[1]*R2[1],m[2]*R2[2],m[3]*R2[3]]; |
> | dotprod(mR2,P); #check that E(mR2)=1 |
> | Em:=dotprod(m,P); #E(m)=1/Rf |
> | sm:=sqrt(P[1]*(m[1]-Em)^2+P[2]*(m[2]-Em)^2+P[3]*(m[3]-Em)^2);#stdev(m) |
> | evalf(%); |
> | (1-Em*ER1)/(sm*s1);#rho(m,R1) |
> | evalf(%); |
> | (1-Em*ER2)/(sm*s2);#rho(m,R2) |
> | evalf(%); |
> | yu:=1/Em+sm/Em*x;#upper part of MV frontier yu=Rf+sigma(m)/E(m)*x |
> | yl:=1/Em-sm/Em*x;#lower part of MV frontier yl=Rf-sigma(m)/E(m)*x |
> | Ex:=ER1*w+(1-w)*ER2;#E(Rp), Rp=wR1+(1-w)R2 |
> | sigma:=sqrt(w^2*s1^2+(1-w)^2*s2^2+2*w*(1-w)*(ER1R2-ER1*ER2));#stdev(Rp) |
> | plot([[x,yu(x),x=0..0.2],[x,yl(x),x=0..0.2],[sigma(w),Ex(w),w=0..1]],color=[red,red,blue]); |
> | msq:=[(12/25)^2, (3/5)^2, (1/2)^2];#calculate m^2 |
> | dotprod(msq,P);#E(m^2) |
> | evalf(%); |
> | moverEmsq:=625/157*[12/25, 3/5, 1/2]; |
> | dotprod(moverEmsq,P);#E(m/E(m^2)) |
> | evalf(%); |
> | mxmoverEmsq:=[m[1]*moverEmsq[1],m[2]*moverEmsq[2],m[3]*moverEmsq[3]]; |
> | dotprod(mxmoverEmsq,P);#check E(m*m/E(m^2))=1 |
> | v1:=[1,0,-6/5]; |
> | u2:=[1,-4,0]; |
> | v1u2:=[1,0,0]; |
> | v1sq := [1, 0, 36/25]; |
> | v2 := u2-dotprod(v1u2, P)*v1/dotprod(v1sq, P); |
> | v1v2 := [144/269, -4*0, -6/5*150/269]; |
> | dotprod(P, v1v2); |
> | v2sq := [20736/72361, 16, 22500/72361]; |
> | Restar := dotprod(P, v1)*v1/dotprod(P, v1sq)+dotprod(P, v2)*v2/dotprod(P, v2sq); |
> | plot3d(-1/4*x2-5/6*x3, x2 = (-10 .. 10), x3 = (-10 .. 10),axes=NORMAL); |
> | [2, 2, 2] = moverEmsq+ga*Restar; |
> | g := (2-300/157)/(7/157); |
> | (2-375/157)/(-61/314); |
> | (2-625/314)/(3/628); |
> | dst:=diff(sigma,w); |
> | dE:=diff(Ex,w); |
> | dEoverdst := proc (w) options operator, arrow; 174*sqrt(4898091*w^2+1269296-4848088*w)/(9796182*w-4848088) end proc; |
> | plot({dEoverdst(w),sqrt(3)/25},w=.5..1); |
> | solve(dEoverdst(w)=sqrt(3)/25); |
> | evalf(%); |
> |