myConstrains[constrains_] := 
 Block[{num, s, mycon}, num = Length[constrains];
  mycon = 
   Table[If[Head[constrains[[i]]] === LessEqual, 
     constrains[[i]][[1]] + s[i] == constrains[[i]][[2]], 
     If[Head[constrains[[i]]] === Equal, 
      constrains[[i]][[1]] + R[i] == constrains[[i]][[2]], 
      constrains[[i]][[1]] - s[i] + R[i] == 
       constrains[[i]][[2]]]], {i, num}];
  mycon /. x__[y_] :> ToExpression[ToString[x] <> ToString[y]]]

LPstd[obj_, mycons_] := 
 Block[{myvar, basicvars, nonbasicvars, AB, BM, c, A, B, rowindex1, 
   rowindex2}, myvar = Variables[mycons[[All, 1]]];
  myvar = 
   Flatten[{Sort@
      Select[myvar, Not@StringStartsQ[SymbolName[#], "s" | "R"] &], 
     Sort@Select[myvar, StringStartsQ[SymbolName[#], "s"] &], 
     Sort@Select[myvar, StringStartsQ[SymbolName[#], "R"] &]}];
  AB = {Coefficient[mycons[[All, 1]], #] & /@ 
     Select[myvar, Not@StringStartsQ[SymbolName[#], "s" | "R"] &], 
    Coefficient[mycons[[All, 1]], #] & /@ 
     Select[myvar, StringStartsQ[SymbolName[#], "s"] &], 
    Coefficient[mycons[[All, 1]], #] & /@ 
     Select[myvar, StringStartsQ[SymbolName[#], "R"] &]};
  AB = Transpose@Flatten[AB, 1];
  AB = Join[{Flatten[
      Table[Select[myvar, 
        StringStartsQ[SymbolName[#], i] &], {i, {"x", "s", "R"}}]]}, 
    AB];
  c = Flatten[{Table[
      Coefficient[obj, #] & /@ 
       Select[myvar, StringStartsQ[SymbolName[#], i] &], {i, {"x", 
        "s"}}], Table[-BM, {Length[
        Select[myvar, StringStartsQ[SymbolName[#], "R"] &]]}], 0}];
  basicvars = 
   Select[AB[[All, #]] & /@ (Flatten /@ 
        Map[Position[myvar, #] &, 
         Tuples[myvar, {Length[mycons]}], {2}]), #[[2 ;; -1]] === 
       IdentityMatrix[Length[mycons]] &][[1, 1]];
  nonbasicvars = 
   Flatten[Table[
     Select[Complement[myvar, basicvars], 
      StringStartsQ[SymbolName[#], i] &], {i, {"x", "s", "R"}}]];
  c = Flatten@{c[[Flatten[Position[myvar, #] & /@ nonbasicvars]]], 
     c[[Flatten[Position[myvar, #] & /@ basicvars]]], 0};
  AB = Join[{Flatten@{nonbasicvars, basicvars}}, 
    Flatten /@ 
     Transpose@{AB[[2 ;; -1]][[All, 
         Sort[Flatten[Position[myvar, #] & /@ nonbasicvars]]]], 
       IdentityMatrix[Length[mycons]]}];
  rowindex1 = c[[-Length[basicvars] - 1 ;; -2]];
  rowindex2 = Range[Length[mycons]];
  Flatten[{c - 
       Total[rowindex1[[#]]*
           Flatten[{AB[[2 ;; -1]][[rowindex2[[#]]]], 
             mycons[[rowindex2[[#]], 2]]}] & /@ 
         Range[Length@rowindex1]], 
      AB[[2 ;; -1]][[All, 
        Sort[Flatten[Position[myvar, #] & /@ nonbasicvars]]]], 
      IdentityMatrix[Length[mycons]], mycons[[All, 2]], 
      Flatten@{Sort@
         Select[nonbasicvars, StringStartsQ[SymbolName[#], "x"] &], 
        Sort@Select[nonbasicvars, 
          StringStartsQ[SymbolName[#], "s"] &], 
        Sort@Select[nonbasicvars, 
          StringStartsQ[SymbolName[#], "R"] &]}, basicvars} & /@ 
    Range[2], 1]]

simplexmethod[c_, A_, B_, b_, xnb_, xb_, cp_, Ap_, Bp_, bp_, xnbp_, 
  xbp_] := 
 Block[{cB, nA, nB, mysymbol, pivotCol, AB, leaveVarpos, Bnew, xbnew, 
   xnbnew, Anew, cnew, znew, z, BM, ABp},
   (*初始化變數*)
   nA = Length[xnb];
   nB = Length[xb];
   
  (*建立符號變數*)
  mysymbol = Flatten@{xnbp, xbp};

  (*找出進入變數（最負的檢驗數）*)
  pivotCol = PositionSmallest[-c[[1 ;; -2]] /. BM -> 10^6][[1]];
  
  (*合併A和B矩陣*)
  AB = Flatten /@ Transpose[{A, B}];
  Print@MatrixForm@AB;
  ABp = Flatten /@ Transpose[{Ap, Bp}];
  leaveVarpos = 
   PositionSmallest[
     b/(If[# <= 0, 10^-6, #] & /@ AB[[All, pivotCol]] /. BM -> 10^6) /. 
      x_ /; x <= 0 -> 10^6][[1]];
  
  (*更新基變數和非基變數*)
  xbnew = xb /. xb[[leaveVarpos]] -> mysymbol[[pivotCol]];
  xnbnew = xnb /. mysymbol[[pivotCol]] -> xb[[leaveVarpos]];
  Bnew = 
   ABp[[All, Flatten[Position[mysymbol, #] & /@ Flatten@{xbnew}]]];
  cB = cp[[Flatten[Position[mysymbol, #] & /@ xbnew]]];
  
  (*計算新的檢驗數和目標函數值*)
  znew = Flatten[{cB . Inverse[Bnew] . Ap - cp[[1 ;; nA]], 
      cB . Inverse[Bnew], cB . Inverse[Bnew] . bp - cp[[-1]]}] // 
    FullSimplify;
  
  (*顯示單形法表格*)
  Print[Grid[
    Insert[Flatten /@ 
      Transpose[{Flatten[{z, xbnew}], 
        Insert[Map[Flatten, 
          Transpose[{Inverse[Bnew] . Ap, Inverse[Bnew], 
            Inverse[Bnew] . bp}]], znew, 1]}], 
     Flatten[{Item["BasicVariable", 
        Background -> Lighter[Yellow, 0.9]], mysymbol, "Solution"}], 
     1], Frame -> All, 
    Background -> {{Lighter[Gray, 0.8]}, {Lighter[Yellow, 0.9], 
       Lighter[Pink, 0.9]}}, ItemStyle -> Directive[FontSize -> 16], 
    Spacings -> {Table[5, {Length[mysymbol] + 3}], 2}]];
  
  (*返回更新後的參數*)
  {-Flatten@{cB . Inverse[Bnew] . Ap - cp[[1 ;; nA]], 
       cB . Inverse[Bnew]}, Inverse[Bnew] . Ap, Inverse[Bnew], 
    Inverse[Bnew] . bp, xnbnew, xbnew, cp, Ap, Bp, bp, xnbp, xbp} // 
   FullSimplify]

Nest[simplexmethod @@ # &, 
 LPstd[-15 x1 - 20 x2, 
  myConstrains[{x1 + 2 x2 >= 10, 2 x1 - 3 x2 <= 6, x1 + x2 >= 6}]], 2]

繼續閱讀全文　Mathematica教學：Simplex method in tableau form

• 編輯文章
• 此文章網址
• 轉寄此文章
• 目前有0 意見

You can download searoute from https://github.com/eurostat/searoute .

port = {
 {"Shanghai, China", 31.2304, 121.4737},
 {"Guangzhou, China", 23.1291, 113.2644},
 {"Rotterdam, Netherlands", 51.9225, 4.4792},
 {"Antwerp, Belgium", 51.2194, 4.4025},
 {"Valencia, Spain", 39.4699, -0.3763},
 {"Los Angeles, USA", 33.7405, -118.2719},
 {"New York/New Jersey, USA", 40.6687, -74.0451},
 {"Savannah, USA", 32.0809, -81.0912},
 {"Santos, Brazil", -23.9608, -46.3336},
 {"Buenos Aires, Argentina", -34.6037, -58.3816},
 {"Itajai, Brazil", -26.9101, -48.6703},
 {"Tangier Med, Morocco", 35.8844, -5.5036},
 {"Port Said, Egypt", 31.2653, 32.3019},
 {"Durban, South Africa", -29.8587, 31.0218},
 {"Melbourne, Australia", -37.8136, 144.9631},
 {"Auckland, New Zealand", -36.8485, 174.7633},
 {"Brisbane, Australia", -27.4698, 153.0251},
 {"Singapore, Singapore", 1.3521, 103.8198},
 {"Kaohsiung, Taiwan", 22.6273, 120.3014},
 {"Busan, South Korea", 35.1796, 129.0756},
 {"Hong Kong, China", 22.3193, 114.1694},
 {"Tanjung Pelepas, Malaysia", 1.3673, 103.5386},
 {"Jebel Ali, UAE", 25.015, 55.0633},
 {"Colombo, Sri Lanka", 6.9271, 79.8612},
 {"Felixstowe, UK", 51.9556, 1.3111},
 {"Algeciras, Spain", 36.1408, -5.4562},
 {"Vancouver, Canada", 49.2827, -123.1207},
 {"Port Klang, Malaysia", 3.0089, 101.3678},
 {"Laem Chabang, Thailand", 13.0846, 100.8956},
 {"Tokyo, Japan", 35.6828, 139.7595},
 {"Yokohama, Japan", 35.4437, 139.6380}
};

Export["/Users/yourname/Desktop/searoute/text.csv",
Insert[
Flatten@{#[[1,1]]<>" to "<>#[[2,1]],
Round[#[[1,3]],0.1],Round[#[[1,2]],0.1],
Round[#[[2,3]],0.1],Round[#[[2,2]],0.1]}&/@
Select[
Flatten[Table[{port[[i]],port[[j]]},
{i,Length@port},{j,Length@port}],1],
#[[1]]!=#[[2]]&
],
{"routename","olon","olat","dlon","dlat"},1]
];

RunProcess[{"java","-jar",
"/Users/yourname/Desktop/searoute/searoute.jar",
"/Users/yourname/Desktop/searoute/text.csv",
"-res","5","-panama","1","-suez","1"}]

route=Import["/Users/yourname/Desktop/searoute/out.geojson",
"Data"][[2,2]];

mypath=Flatten@{#[[2,2]],StringSplit["routename"/.#[[3,2]]," to "],
"distKM"/.#[[3,2]]}&/@route;

myrule=Rule@@@Transpose@{Tally[mypath[[All,2]]][[All,1]],
Range[Length@port]};

distMatrix=ParallelTable[If[i==j,Infinity,
ToExpression@Select[mypath/.myrule,
#[[2]]==i&&#[[3]]==j&][[1,-1]]],
{i,Length@port},{j,Length@port}];

tour=FindShortestTour[Range[Length@port],
DistanceFunction->(distMatrix[[#1,#2]]&)];

myline=ParallelTable[
If[i==j,{},Select[mypath/.myrule,#[[2]]==i&&#[[3]]==j&][[1,1]]],
{i,Length@port},{j,Length@port}];

myline=myline[[#[[1]],#[[2]]]]&/@Partition[tour[[2]],2,1];

GeoGraphics[{
Red,PointSize[Large],
Point[GeoPosition[port[[All,{2,3}]]]],
Blue,Thick,myline},
GeoRange->"World",
GeoCenter->GeoPosition[{20,-150}],
GeoRange->Quantity[180,"AngularDegrees"]
]

繼續閱讀全文　Mathematica 教學：Searoute for Traveling Salesman Problem

• 編輯文章
• 此文章網址
• 轉寄此文章
• 目前有0 意見

Concise and elegant!

g=Graph[
{1\[DirectedEdge]2,1\[DirectedEdge]3,2\[DirectedEdge]3,
2\[DirectedEdge]4,3\[DirectedEdge]5,4\[DirectedEdge]3,
4\[DirectedEdge]5,4\[DirectedEdge]6,5\[DirectedEdge]6},
EdgeWeight->{1,12,9,3,5,3,13,15,4},
EdgeLabels->Rule@@@Transpose@{
{1\[DirectedEdge]2,1\[DirectedEdge]3,2\[DirectedEdge]3,
2\[DirectedEdge]4,3\[DirectedEdge]5,4\[DirectedEdge]3,
4\[DirectedEdge]5,4\[DirectedEdge]6,5\[DirectedEdge]6},
{1,12,9,3,5,3,13,15,4}},
VertexStyle->Yellow,VertexSize->0.5,
VertexShapeFunction->{{Yellow,Disk[#1,0.1],Black,Text[#2,#1]}&}]

route={{1,2,1},{1,3,12},{2,3,9},{2,4,3},{3,5,5},{4,3,4},{4,5,13},{5,6,4},{4,6,15}};

myroute=Select[Table[Select[route,#[[2]]==z&],{z,1,6}],#=!={}&]

Table[Set[v[temp], 
SortBy[If[#[[1]]==1,{v[temp],#[[3]]},{v[temp],v[#[[1]]]+#[[3]]}]&/@
Select[Flatten[myroute,1],#[[2]]==temp&],Last][[1,2]]],{temp,2,6}]&/@Range[2]

繼續閱讀全文　Mathematica 教學：Dijkstra Algorithm

• 編輯文章
• 此文章網址
• 轉寄此文章
• 目前有0 意見

A Sustainable Dynamic Optimization Model of Pricing and Advertising in the Presence of Green Innovation Investment

Chung-Yuan Dye and Tsu-Pang Hsieh

Finally, the acceptance from EJOR has arrived. The Julia code for the numerical examples in this paper is provided below. Alternatively, you can also implement it directly in Julia. Please leave a message if you have any problems with the implementation of the code.

julia=StartExternalSession["Julia"]

ExternalEvaluate[julia,"
using JuMP,Ipopt
dp(x)=(140-3*x);
K1(y)=(y/5-y^2/1000);
f1(x)=100*exp(-0.03*x);
f2(x)=2*exp(-0.045*x);
v1=25;
v2=10;
omega=300;
c=16.5;
h=0;
z=0.8;
"];

dynCP11=ExternalEvaluate[julia, "
	function dynCP11(G0,gamma1,gamma2,delta,k1coef,k2coef,e,r,varpi,n)
    m1 = Model(Ipopt.Optimizer) 
    set_silent(m1)
    #set_optimizer_attributes(m1, \"tol\" => 1e-4, \"max_iter\" => 100)
    register(m1, :dp, 1, dp; autodiff = true)
    register(m1, :K1, 1, K1; autodiff = true)
    register(m1, :f1, 1, f1; autodiff = true)
    register(m1, :f2, 1, f2; autodiff = true)
    @variable(m1, p[1:n], start = 30)
    @variable(m1, w[1:n], start = 10)
    @variable(m1, u[1:n], start = 10)

    @NLexpression(m1,newG[i=1:n],(1-delta)^(-1+i)*((1-delta)*G0
	+sum((gamma1*u[j+1]^k1coef+gamma2*w[j+1]^k2coef)/((1-delta)^j) 
    for j = 0:i-1)))
    
    @NLexpression(m1,newCE[i=1:n],f1(w[i])+f2(w[i])*dp(p[i])*K1(newG[i])-varpi)

    @NLobjective(
	m1,
	Max,
	sum(z^(i-1)*((p[i]-c-h/2)*dp(p[i])*K1(newG[i])
    -v1*u[i]^2/2-v2*w[i]^2/2-e*newCE[i])
                for i = 1:n))
	for i = 1:n
        @constraint(m1, p[i] >= c + h / 2)
	end
	for i = 1:n
        @constraint(m1, p[i] <= 140 / 3)
	end
	for i = 1:n
        @constraint(m1, w[i] >= 0)
	end
	for i = 1:n
        @constraint(m1, u[i] >= 0)
	end
	for i = 1:n
        @NLconstraint(m1, newCE[i] >= 0)
	end

	JuMP.optimize!(m1)

	p_val = JuMP.value.(p)
	w_val = JuMP.value.(w)
	u_val = JuMP.value.(u)

	G_val=(1-delta)*G0+gamma1*u_val[1]^k1coef+gamma2*w_val[1]^k2coef

	CE=f1(w_val[1])+f2(w_val[1])*dp(p_val[1])*K1(G_val)

	profit=(p_val[1]-c-h/2)*dp(p_val[1])*K1(G_val)
    -v1*u_val[1]^2/2-v2*w_val[1]^2/2-e*(CE-varpi)

	return[p_val,w_val,u_val,G_val,CE,profit,objective_value(m1)]

end"];

dynCP22=ExternalEvaluate[julia,"
	function dynCP22(G0,gamma1,gamma2,delta,k1coef,k2coef,e,r,varpi,n)
    m1 = Model(Ipopt.Optimizer) 
    set_silent(m1)
    #set_optimizer_attributes(m1, \"tol\" => 1e-4, \"max_iter\" => 100)
    register(m1, :dp, 1, dp; autodiff = true)
    register(m1, :K1, 1, K1; autodiff = true)
    register(m1, :f1, 1, f1; autodiff = true)
    register(m1, :f2, 1, f2; autodiff = true)
    @variable(m1, p[1:n], start = 30)
    @variable(m1, w[1:n], start = 10)
    @variable(m1, u[1:n], start = 10)

    @NLexpression(m1,newG[i=1:n],(1-delta)^(-1+i)*((1-delta)*G0
	+sum((gamma1*u[j+1]^k1coef+gamma2*w[j+1]^k2coef)/((1-delta)^j) 
    for j = 0:i-1)))
    
    @NLexpression(m1,newCE[i=1:n],f1(w[i])+f2(w[i])*dp(p[i])*K1(newG[i])-varpi)

    @NLobjective(
	m1,
	Max,
	sum(z^(i-1)*((p[i]-c-h/2)*dp(p[i])*K1(newG[i])
    -v1*u[i]^2/2-v2*w[i]^2/2-r*newCE[i])
                for i = 1:n))
	for i = 1:n
        @constraint(m1, p[i] >= c + h / 2)
	end
	for i = 1:n
        @constraint(m1, p[i] <= 140 / 3)
	end
	for i = 1:n
        @constraint(m1, w[i] >= 0)
	end
	for i = 1:n
        @constraint(m1, u[i] >= 0)
	end
	for i = 1:n
        @NLconstraint(m1, newCE[i] <= 0)
	end

	JuMP.optimize!(m1)

	p_val = JuMP.value.(p)
	w_val = JuMP.value.(w)
	u_val = JuMP.value.(u)

	G_val=(1-delta)*G0+gamma1*u_val[1]^k1coef+gamma2*w_val[1]^k2coef

	CE=f1(w_val[1])+f2(w_val[1])*dp(p_val[1])*K1(G_val)

	profit=(p_val[1]-c-h/2)*dp(p_val[1])*K1(G_val)
    -v1*u_val[1]^2/2-v2*w_val[1]^2/2-r*(CE-varpi)

	return[p_val,w_val,u_val,G_val,CE,profit,objective_value(m1)]

end"];

mysol13[g0_,\[Gamma]1_,\[Gamma]2_,\[Delta]_,k1_,k2_,e_,r_,cap_]:=
Block[{n=40,numiter=100,iter=0,ans1,ans2},
	Print[{Dynamic[iter],Dynamic[ans2]}];
	NestList[(
	iter=iter+1;
	ans1=
	SortBy[{
		dynCP11[#[[4]],\[Gamma]1,\[Gamma]2,\[Delta],k1,k2,e,r,cap,n],
		dynCP22[#[[4]],\[Gamma]1,\[Gamma]2,\[Delta],k1,k2,e,r,cap,n]},
          Last][[-1]];
	ans2={ans1[[1,1]],ans1[[2,1]],ans1[[3,1]],ans1[[4]],
          ans1[[5]],#[[-1]]+0.8^(iter-1)*ans1[[6]]})&,
	{0,0,0,g0,0,0},numiter][[2;;-1]]
]

myexample=mysol13[#[[1]],0.5,0.8,0.2,0.5,0.8,8,2,#[[2]]]&/@
Flatten[Table[{i,j},{j,{250,300,400}},{i,{10,40}}],1]

ListLinePlot[myexample[[All,All,#]],PlotRange->All,
Frame->True]&/@{1,2,3,4,5}//Column

繼續閱讀全文　Finally, the first acceptance of 2024 comes.

• 編輯文章
• 此文章網址
• 轉寄此文章
• 目前有0 意見

test=ImageData[Rasterize["50",RasterSize->50,ImageSize->200]]

data=ParallelTable[{i,j,If[Mean@test[[i,j]]>0.95,0,1]},
{i,Dimensions[test][[1]]},{j,Dimensions[test][[2]]}]

picdata=Table[ParallelMap[If[#[[-1]]==0,"",
ImageResize[Import[RandomChoice[FileNames["*.jpeg",
"/Users/chungyuandye/Pictures/照片圖庫V4.photoslibrary/originals/",Infinity]]],{100,100}]]&,
data[[i]]],{i,54}]

GraphicsGrid[picdata]

繼續閱讀全文　Mathematica 教學：Text arrangement

• 編輯文章
• 此文章網址
• 轉寄此文章
• 目前有0 意見

session=StartExternalSession["Julia"]
test =
  "function testpi(size)
    sample = rand(size, 2)
    dist = broadcast(hypot, sample[:,1], sample[:,2])
    return sample, dist
end";
pitest=ExternalEvaluate[session,test]
size=100000
data=pitest[size];
ListPlot[GatherBy[Transpose[data[[1]]],Norm[#]<1&],Frame ->True,AspectRatio -> 1]
Length@GatherBy[Transpose[data[[1]]], Norm[#]<1&][[1]]/size//N




\documentclass[tikz]{standalone}
\usepackage{tikz}
\usepackage{ifthen}
\begin{document}
  \begin{tikzpicture}
	\draw[help lines, color=gray!30, dashed] (-1.2,1.5) grid (-1.2,1.5);
	\draw[->] (-1.2,0)--(1.1,0) node[right]{$x$};
	\draw[->] (0,-1.2)--(0,1.5) node[above]{$y$};
	\draw[domain=0:1, smooth, variable=\x, black] plot ({\x}, {sqrt(1-\x*\x)});
  \pgfmathsetseed{42}
    \foreach \i in {1,...,10000}
    {
      \pgfmathsetmacro{\x}{abs(rand)}
      \pgfmathsetmacro{\y}{abs(rand)}
      \pgfmathsetmacro{\r}{\x*\x+\y*\y}
       \pgfmathsetmacro{\col}{ifthenelse(\r<1,"blue!90","red!90")}
      \fill[\col] (\x,\y) circle (0.1pt)
      ;
    }
  \end{tikzpicture}
\end{document}

繼續閱讀全文　Mathematica 教學：Call Julia from Mathematica

• 編輯文章
• 此文章網址
• 轉寄此文章
• 目前有3 意見

 

Policy function iterations方法是通過計算出最佳化問題的每一步所實現的最佳結果，並將這些結果用於計算下一步的最佳結果。此過程在一定數量的迭代後會得到一個最佳策略，並且這個策略也是最終的最佳解。因此它對於解決複雜的最佳化問題非常有用。它可以被用於解決各種經濟學、工程學和其他領域的最佳化問題。

ClearAll["Global`*"]
\[Alpha]=0.3;
\[Beta]=0.95;
\[Delta]=0.1;
\[Sigma]=1.5;

size=99;

K=Range[0.1,5,4.9/(size)];
V=Table[0,{size}];


astHoward=With[{n=size},
obj=ParallelTable[If[i^\[Alpha]+(1-\[Delta])*i-j<0,-10^9,
((i^\[Alpha]+(1-\[Delta])*i-j)^(1-\[Sigma])-1)/(1-\[Sigma])],
{i,0.1,5,4.9/n},{j,0.1,5,4.9/n}];
loops=0;
\[Epsilon]=99999;

Print[Grid[{
{"loop number ",
Dynamic[loops],
" with distance ",
Dynamic[NumberForm[\[Epsilon],
ExponentFunction->(If[-10<#<10,Null,#]&)]]}}]];

NestWhile[
NestWhile[
Transpose@
Table[{obj[[i,#[[2,i]]]]+\[Beta]*#[[1,#[[2,i]]]],#[[2,i]]},{i,Length[#[[1]]]}]&,
Transpose[SortBy[#,First][[-1]]&/@
Table[{obj[[i,j]]+\[Beta]*#[[1,j]],j},{i,n+1},{j,n+1}]//Parallelize],
Norm[#1[[1]]-#2[[1]]]>10^-6&,2]&,
{Table[0,n+1],Table[0,n+1]},
(loops++;\[Epsilon]=Norm[#1[[1]]-#2[[1]]];
Norm[#1[[1]]-#2[[1]]]>10^-6)&,2]
]


astHowardMatrix=With[{n=size},
obj=ParallelTable[If[i^\[Alpha]+(1-\[Delta])*i-j<
0,-10^9,((i^\[Alpha]+(1-\[Delta])*i-j)^(1-\[Sigma])-1)/(1-\[Sigma])],
{i,0.1,5,4.9/n},{j,0.1,5,4.9/n}];
loops=0;
\[Epsilon]=99999;

Id=IdentityMatrix[size+1];

Print[Grid[{
{"loop number ",
Dynamic[loops],
" with distance ",
Dynamic[NumberForm[\[Epsilon],
ExponentFunction->(If[-10<#<10,Null,#]&)]]}}]];


NestWhile[
(
myValue=Transpose[SortBy[#,First][[-1]]&/@
Table[{obj[[i,j]]+\[Beta]*#[[1,j]],j},{i,n+1},{j,n+1}]//Parallelize];
{Inverse[(Id-\[Beta]*Normal[SparseArray[
Rule[#,1]&/@Transpose[
{Range[size+1],myValue[[2]]}],size+1]])].(obj[[#[[1]],#[[2]]]]&/@
Transpose[{Range[size+1],myValue[[2]]}]),myValue[[2]]}
)&,
{Table[0,n+1],Table[0,n+1]},
(loops++;\[Epsilon]=Norm[#1[[1]]-#2[[1]]];
Norm[#1[[1]]-#2[[1]]]>10^-6)&,2]
]


ListLinePlot[astHowardMatrix[[1]],PlotRange->All,AspectRatio->1,Frame->True]
ListLinePlot[astHowardMatrix[[2]],PlotRange->All,AspectRatio->1,Frame->True]

繼續閱讀全文　Mathematica 教學：Policy function iterations for deterministic cake eating problem

• 編輯文章
• 此文章網址
• 轉寄此文章
• 目前有0 意見

Value function iterations方法是通過計算出最佳化問題的每一步所實現的最佳結果，並將這些結果用於計算下一步的最佳結果。此過程在一定數量的迭代後會得到一個最佳策略，並且這個策略也是最終的最佳解。因此它對於解決複雜的最佳化問題非常有用。它可以被用於解決各種經濟學、工程學和其他領域的最佳化問題。

ClearAll[\[Beta],\[Gamma],xGrid,u,v0];

(*Parameters*)
\[Beta]=0.96;
\[Gamma]=1.5;
xMin=10^-3;
xMax=2.5;
nGrid=120;
xGrid=Subdivide[xMin,xMax,nGrid-1];

(*Utilityfunction*)
u[c_]:=If[\[Gamma]==1,Log[c],(c^(1-\[Gamma]))/(1-\[Gamma])];

(*Initialguess*)
v0=Table[0,{x,xGrid}];

(* Bellman Operator T *)
myT[v_List]:=Block[{myf,x},
myf[x_]=Interpolation[Transpose@{xGrid,v},InterpolationOrder->1][x];
Max/@Quiet@Table[
(c^(1-\[Gamma]))/(1-\[Gamma])+\[Beta]*myf[x-c],
{x,xGrid},{c,Select[xGrid,#<=x&]}
]];

(* Value Function Iteration using NestWhile *)
myvfi=NestWhileList[myT,v0,Norm[#1-#2]>=0.00001&,2];

(* Plotting the results *)
ListLinePlot[Transpose[{xGrid,#}]&/@myvfi[[2;;-1]]]

繼續閱讀全文　Mathematica 教學：Value function iterations for deterministic cake eating problem

• 編輯文章
• 此文章網址
• 轉寄此文章
• 目前有0 意見

splitdata[datatemp___]:=Block[{eqn,temp,x,y,group,vertex,ans},
temp=SortBy[datatemp,Left][[{1,-1}]];
eqn[x_,y_]=y-temp[[1,2]]-
(temp[[1,2]]-temp[[2,2]])/(temp[[1,1]]-temp[[2,1]])(x-temp[[1,1]]);
group={Select[datatemp,eqn[#[[1]],#[[2]]]> 0&],
Select[datatemp,eqn[#[[1]],#[[2]]]< 0&]};

vertex={SortBy[group[[1]],Abs@eqn[#[[1]],#[[2]]]&][[-1]],
SortBy[group[[2]],Abs@eqn[#[[1]],#[[2]]]&][[-1]]};

ans={{{temp[[1]],vertex[[1]]},
Select[group[[1]],RegionMember[
Polygon[Insert[temp,vertex[[1]],-1]],#]==False&&#[[1]]< vertex[[1,1]]&]},
{{temp[[2]],vertex[[1]]},
Select[group[[1]],RegionMember[
Polygon[Insert[temp,vertex[[1]],-1]],#]==False&&#[[1]]> vertex[[1,1]]&]},
{{temp[[1]],vertex[[2]]},
Select[group[[2]],RegionMember[
Polygon[Insert[temp,vertex[[2]],-1]],#]==False&&#[[1]]< vertex[[2,1]]&]},
{{temp[[2]],vertex[[2]]},
Select[group[[2]],RegionMember[
Polygon[Insert[temp,vertex[[2]],-1]],#]==False&&#[[1]]> vertex[[2,1]]&]}
};

If[Length[#[[2]]]==0,{#[[1]],{#[[1,1]]}},{#[[1]],#[[2]]}]&/@ans
];

convexhull[{p1__,p2__},datatemp___]:=Block[{eqn,temp,x,y,group,vertex,ans},
temp=SortBy[{p1,p2},Left];
If[Length[datatemp]==0,{{temp,{}},{temp,{}}},
If[Length[datatemp]==1,{{{temp[[1]],datatemp[[1]]},{}},
{{temp[[2]],datatemp[[1]]},{}}},

eqn[x_,y_]=y-temp[[1,2]]-(temp[[1,2]]-temp[[2,2]])/
(temp[[1,1]]-temp[[2,1]])(x-temp[[1,1]]);

vertex=SortBy[datatemp,Abs@eqn[#[[1]],#[[2]]]&][[-1]];

group={Select[datatemp,#[[1]]< vertex[[1]]&],
Select[datatemp,#[[1]] > vertex[[1]]&]};

ans={{{temp[[1]],vertex},
Select[group[[1]],RegionMember[Polygon[Insert[temp,vertex,-1]],#]==
False&&#[[1]]< vertex[[1]]&]},
{{vertex,temp[[2]]},
Select[group[[2]],RegionMember[Polygon[Insert[temp,vertex,-1]],#]==
False&&#[[1]]> vertex[[1]]&]}}
]]
];

data=RandomReal[NormalDistribution[0,2],{250,2}];
Needs["ComputationalGeometry`"]

test=splitdata[data];
xx=Nest[Flatten[convexhull@@@#,1]&,test,3];
myline=Line/@Tally[xx[[All,1]]][[All,1]];
mypts=Point@Tally[xx[[All,1]]][[All,1,1]];

ListPlot[data,AspectRatio->1,
Epilog->{Red,Thickness[0.015],myline,Blue,PointSize[0.05],mypts}]
 

繼續閱讀全文　Mathematica 教學：Convex Hull

• 編輯文章
• 此文章網址
• 轉寄此文章
• 目前有0 意見

  

論文統計分析免煩惱，平均數比較、迴歸分析、信度、效度 CFA、結構方程是模型 SEM，一鍵搞定！ 

範例：本研究採量化研究，以XXX為研究對象，發放問卷XXX份，回收有效問卷XXX份，並以三山出版社所提供的 ThesisNoWorries 論文統計分析免煩惱 統計套件進行描述性統計、信效度分析、平均數比較、迴歸和相關分析、及結構方程式模型分析。

程式碼：https://github.com/chungyuandye/ThesisNoWorries

繼續閱讀全文　ThesisNoWorries Beta 0.3

• 編輯文章
• 此文章網址
• 轉寄此文章
• 目前有0 意見

mytwoaxistikz[data__,{xmin_,xmax_,ymin_,ymax_},
mycolor__,{xlabel_,ylabel_},{rticks__,decimalr_,lticks__,decimall_},
linelabel__,linepos__,plotaspectratio_:1]:=Block[{temp=data,axis1,axis2,textpos,line},
axis1="
\\begin{axis}[
    yticklabel style={
        /pgf/number format/.cd,
            fixed,
            fixed zerofill,
            precision="<>ToString[decimall]<>",
    },
 scale only axis, 
 axis y line*=left, 
 xtick pos=left, 
 xscale="<>ToString[N@plotaspectratio]<>", 
 xmin="<>ToString[xmin]<>",xmax="<>ToString[xmax]<>",ymin="<>ToString[ymin]<>
 ",ymax="<>ToString[ymax]<>",
 xlabel="<>xlabel<>",ylabel="<>ylabel<>",
 ytick="<>ToString[lticks]<>"]\n";

textpos=StringJoin[Table["
 \\node[] at (axis cs:"<>ToString[linepos[[z,1]]]<>","<>ToString[linepos[[z,2]]]<>
 ") {"<>linelabel[[z]]<>"};",{z,Length@temp}]];

line=StringJoin[Table["
\\addplot[color="<>If[ListQ[mycolor[[z]]],
ToString[mycolor[[z,1]]]<>","<>ToString[mycolor[[z,2]]],
ToString[mycolor[[z]]]]<>",mark=none,line width=0.5mm] coordinates{\n"
 <>StringJoin["("<>ToString[#[[1]]]<>","<>ToString[#[[2]]]<>")\n"&/@temp[[z]]]<>
 "};\n\\addlegendentry{$G_{0}= "<>ToString[G0[[z]]]<>"$}",{z,Length@temp}]]
<>"\n\\end{axis}\n";

axis2="\n\\begin{axis}[ 
    yticklabel style={
        /pgf/number format/.cd,
            fixed,
            fixed zerofill,
            precision="<>ToString[decimalr]<>",
    },
 scale only axis, 
 axis y line*=right, 
 xtick pos=left, 
 xscale="<>ToString[N@plotaspectratio]<>", 
 xmin="<>ToString[xmin]<>",xmax="<>ToString[xmax]<>",ymin="<>ToString[ymin]<>
 ",ymax="<>ToString[ymax]<>",
 ytick="<>ToString[rticks]<>"]\n";

"\\begin{tikzpicture}"<>axis1<>textpos<>line<>axis2<>"\n\\end{axis}\n\\end{tikzpicture}"
]

繼續閱讀全文　Exporting graphics to TikZ format with Mathematica, Version 2

• 編輯文章
• 此文章網址
• 轉寄此文章
• 目前有1 意見

安裝 Ijupyter

python3 -m pip install --upgrade pip
python3 -m pip install jupyter

設定遠端連線 jupyter server，在python shell中執行

from notebook.auth import passwd; passwd()

生成修改設定檔

jupyter notebook --generate-config
vi ~/.jupyter/jupyter_notebook_config.py
c.NotebookApp.ip='*' #*表示所有IP，這裡設定ip為都可訪問  
c.NotebookApp.password = '填寫剛剛生成的密文'
c.NotebookApp.open_browser = False # 禁止notebook啟動時自動開啟瀏覽器
c.NotebookApp.port =9999 #指定訪問的埠，預設是8888。 

下載 configure-jupyter.wls

https://github.com/WolframResearch/WolframLanguageForJupyter

執行 wolframscript 安裝 WolframLanguageForJupyter

wolframscript configure-jupyter.wls add

啟動 Jupyter notebook

繼續閱讀全文　Mathematica 教學：使用 Jupyter Notebook 做為 Mathematica 操作環境

• 編輯文章
• 此文章網址
• 轉寄此文章
• 目前有0 意見

xx=Import["http://www.fishbase.org/summary/Decapterus-maruadsi.html","Table"];
StringJoin[#<>""&/@Cases[xx,{"Fisheries:",___},Infinity][[1,1;;-2]]]

繼續閱讀全文　Mathematica 教學：抓網頁特定資料

• 編輯文章
• 此文章網址
• 轉寄此文章
• 目前有0 意見

1. MI值調整
2. 項目分析報表修訂
3. 三構面中介模型報表修訂
4. 新增兩構面因果關係路徑模式

繼續閱讀全文　 ThesisNoWorryMed Beta 0.2

• 編輯文章
• 此文章網址
• 轉寄此文章
• 目前有0 意見

mydate[date_]:=Block[{temp},
temp=Import[
"http://www.twse.com.tw/exchangeReport/MI_INDEX.html?date="<>ToString@date<>"&type=ALLBUT0999",
"JSON"];
{("data1"/.temp)[[All,{1,2,4,5}]],("data5"/.temp)[[All,
{1,2,3,4,5,6,7,8,9,11}]]}]

mydate[20170719][[1]] // TableForm
mydate[20170719][[2]] // TableForm

test=Quiet@{
{ToExpression@StringJoin@Characters[ToString@#][[1;;4]],
ToExpression@StringJoin@Characters[ToString@#][[5;;6]],
ToExpression@StringJoin@Characters[ToString@#][[7;;8]]},
mydate[#][[1,4]]}&/@Map[ToExpression@
StringJoin[ToString/@{#[[1]],
If[#[[2]]<10>ToString[#[[2]]],ToString[#[[2]]]],
If[#[[3]]<10>ToString[#[[3]]],
ToString[#[[3]]]]}]&,
DayRange["2017/05/1","2017/07/21"][[All,1]]]//Parallelize;

DateListPlot[{#[[1]],
ToExpression@StringReplace[#[[2,2]],","->""]}&/@
Select[test,Length@#[[2]]==4&],Filling->Bottom]


Block[{tempdata},
tempdata={#[[1]],ToExpression@StringReplace[#[[2,2]],","->""]}&/@
Select[test,Length@#[[2]]==4&];
DateListPlot[tempdata,Filling->Bottom,
Epilog->{Red,PointSize[0.025],Point[{tempdata[[#[[1]]]][[1]],#[[2]]}]&/@
FindPeaks[tempdata[[All,2]]]}]]
 

繼續閱讀全文　Mathematica 教學：抓取證交所股價資料

• 編輯文章
• 此文章網址
• 轉寄此文章
• 目前有0 意見

Needs["RLink`"]

(* For Mac and MMA 10*)

SetEnvironment[ "DYLD_LIBRARY_PATH" ->"/Library/Frameworks/R.framework/Resources/lib"]
InstallR["RHomeLocation" -> "/Library/Frameworks/R.framework/Resources", "RVersion" -> 3];
REvaluate["R.version.string"]

(* Load packages*)

REvaluate["library(psych);"]
REvaluate["library(nortest);"]
REvaluate["library(GenABEL);"]

(* 因素分析，主成分法萃取，變異及大化法轉軸 *)

myfac=Map[REvaluate["{xx<-principal(iris[,1:4],nfactors="<>ToString[#]<>",
rotate=\"varimax\");xx$loadings}"]&,{2,3,4}];
myfac

(* 輸出因素負荷矩陣 *)

Reverse[SortBy[myfac[[1,1]],{#[[1]],#[[2]]}&]]//TableForm
TableForm/@myfac[[All,1]]


(* 因素分析資料*) 

fadata = {
{1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 3, 3, 3, 3, 3, 4, 5, 6}, 
{1, 2, 1, 1, 1, 1, 2, 1, 2, 1, 3, 4, 3, 3, 3, 4, 6, 5}, 
{3, 3, 3, 3, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 5, 4, 6}, 
{3, 3, 4, 3, 3, 1, 1, 2, 1, 1, 1, 1, 2, 1, 1, 5, 6, 4}, 
{1, 1, 1, 1, 1, 3, 3, 3, 3, 3, 1, 1, 1, 1, 1, 6, 4, 5}, 
{1, 1, 1, 2, 1, 3, 3, 3, 4, 3, 1, 1, 1, 2, 1, 6, 5, 4}};

(* 定義資料*)

RSet["factor", Transpose@fadata];

(* 在R中計算相關係數矩陣，傳回MMA小數後四捨五入 *)

Round[REvaluate["cor(factor)"], 0.0001] // MatrixForm

(* 因素分析 *)

REvaluate["factanal(factor, factors = 3)$loadings"]

(* 輸出因素負荷矩陣 *)

TableForm[#[[1]], TableHeadings -> {None, #[[-1, -1, -1, -1]]}] & /@ 
{REvaluate["factanal(factor, factors = 3)$loadings"]} // Column
 
 
(* 迴歸分析 *)

REvaluate["{
data(iris)
reg <- lm( Sepal.Length ~ Species, data=iris )
summary.text <- capture.output(print( summary(reg)) )
}"]

StringJoin@Riffle[#, "\n"] &@REvaluate["{
data(iris)
reg <- lm( Sepal.Length ~ Species, data=iris )
summary.text <- capture.output(print( summary(reg)) )
}"]
 
(* 產生常態隨機變數 *)

mydata=REvaluate["rnorm(100,mean=0,sd=1)"]&/@Range[100]; 

(* MMA次數分配 *)

\[ScriptCapitalD] = HistogramDistribution[mydata[[1]]]  

(* MMA PDF & CDF *)

GraphicsRow[DiscretePlot[#[\[ScriptCapitalD], x], {x, -4, 4, .01}, 
PlotLabel -> #] & /@ {PDF, CDF}]

(* shapiro test, ks.test *)

shapiro[data_] := Block[{temp},
  RSet["temp", data];
  REvaluate[
 "{xx=shapiro.test(temp); yy=ks.test(temp,\"pnorm\",0,1); zz=lillie.test(temp); 
 c(xx$statistic,xx$p.value,yy$statistic,yy$p.value,zz$statistic,zz$p.value)}"][[1]] ]
 
RSet["rn1", RandomReal[{0, 10}, 1000]];
rn2 = RandomReal[NormalDistribution[0, 1], 1000];
REvaluate["rntransform(rn1)"] // Short
Histogram[#, Frame -> True] & /@ {rn2, REvaluate["rntransform(rn1)"]}
ShapiroWilkTest[#, {"TestStatistic", "PValue"}] & /@ {rn2, REvaluate["rntransform(rn1)"]} 


(* 利用RCurl擷取網路資料 *)

gosspost[n_]:=REvaluate["
{library(RCurl);
curl<-getCurlHandle();
curlSetOpt(cookie=\"over18=1\",followlocation=TRUE,curl=curl);
url<-\"https://www.ptt.cc/bbs/Gossiping/index" <> ToString[n] <> 
".html\";
html<-getURL(url,curl=curl)
}"][[1]];

Gossiping = 
  Cases[ImportString[
StringReplace[gosspost[#], {"\n" -> "", "\t" -> ""}], 
"XMLObject"],
  XMLElement[
 "div", {"class" -> "r-ent"}, {XMLElement[
"div", {"class" -> "nrec"}, {XMLElement[
  "span", {"class" -> "hl f2"}, {_}]}], 
  XMLElement["div", {"class" -> "mark"}, {}], 
  XMLElement[
"div", {"class" -> "title"}, {XMLElement[
  "a", {"shape" -> "rect", 
"href" -> postlink_}, {posttile_}]}], 
  XMLElement[
"div", {"class" -> "meta"}, {XMLElement[
  "div", {"class" -> "date"}, {_}], 
 XMLElement[
  "div", {"class" -> "author"}, {author_}]}]}] :> {postlink, 
 posttile, author},
  Infinity] & /@ Range[20861, 20862, 1];
  
Flatten[Gossiping, 1] // TableForm  

繼續閱讀全文　Use R within Mathematica

• 編輯文章
• 此文章網址
• 轉寄此文章
• 目前有0 意見

1+2Sin[2x]+3Sin[3x]+4Sin[4x]/.{Plus->List,Sin[y_]->1}

繼續閱讀全文　How to get the coefficient of polynomial function

• 編輯文章
• 此文章網址
• 轉寄此文章
• 目前有0 意見

str=Reverse@Rasterize[Style[\[Pi],FontFamily->"Kokonor",Bold,Red,FontSize->20], 
RasterSize->20,ImageSize->{Automatic,25}][[1,1]][[6;;-12]];

mytext=Map[If[Total@#==255*3,"",0]&,str,{2}];

mypos=Position[Map[If[Total@#==255*3,1,0]&,str,{2}],0];

mypi=Drop[Characters[ToString@N[Pi,Length@mypos]],{2}];

Set[mytext[[mypos[[#,1]],mypos[[#,2]]]],Style[mypi[[#]],
Lighter@RandomColor[]]]&/@Range[Length@mypi];

Grid[mytext,Spacings->{0.2,0.2},Background->Black,
ItemStyle->{Directive[{RandomReal[255],RandomReal[255],
RandomReal[255]},18,Italic],None}]

繼續閱讀全文　Mathematica 教學：Pi Day

• 編輯文章
• 此文章網址
• 轉寄此文章
• 目前有0 意見

mypic=your pic;
ParametricPlot3D[{v*Cos[u],v*Sin[u],0.0035v},{u,0,2Pi},{v,0,1},
Mesh->0,BoundaryStyle->None,PlotStyle->Texture[mypic],
ViewPoint->{0,0,1},Boxed->False,Axes->None,PlotPoints->50,
ImagePadding->None,
Lighting->{{"Directional",White,{{0,3,2},{0,3,0}}}}]

ParametricPlot3D[{v*Cos[u],v*Sin[u],0.0035v},{u,0,2Pi},{v,0,1},
Mesh->0,BoundaryStyle->None,PlotStyle->Texture[ImageRotate[mypic,Pi]],
ViewPoint->{0,0,1},Boxed->False,Axes->None,PlotPoints->50,
ImagePadding->None,
Lighting->{{"Directional",White,{{0,3,2},{0,3,0}}}}]

繼續閱讀全文　Mathematica 教學：Generate a Tiny Planet Photo

• 編輯文章
• 此文章網址
• 轉寄此文章
• 目前有0 意見

(* 輸入資料 *)
data[i_]:=Import["https://www.ptt.cc/bbs/Mathematica/index"<>
ToString@i<>".html",{"XMLObject"}];

(* 讀取文章欄位資料 *)
mydata=Flatten[Flatten/@
Transpose@
{
Cases[data[#],XMLElement["div",{"class"->"date"},{date_}]:>date,
Infinity],
Cases[data[#],XMLElement["div",{"class"->"author"},{author_}]:>author,
Infinity],
Cases[data[#],XMLElement["a",{"shape"->"rect","href"->link_},{title_}]:>{title,link},
Infinity]
}&/@Range[34],1]

(* 分割年度 *)
mydate=Split[Quiet@Flatten@DateList[StringReplace[#,""->""]][[{2,3}]]&/@
mydata[[All,1]],#1[[1]]<=#2[[1]]&];

(* 重新定義日期  *)
mydata[[All,1]]=Flatten[MapThread[
Flatten/@Thread[{#1,#2}]&,{Range[2011,2016],mydate}],1];

(* 文章直方圖 *)
DateHistogram[mydata[[All,1]],#]&/@{"Month","Year"}

DateHistogram[Select[mydata,#[[2]]=="chungyuandye"&][[All,1]],#]&/@{"Month","Year"}

繼續閱讀全文　Mathematica 教學：讀取PTT資料

• 編輯文章
• 此文章網址
• 轉寄此文章
• 目前有0 意見