Instead of using the linear model y = a + b*x, use y = b*x. If you're using Statistics:-LinearFit, then the command is

Statistics:-LinearFit(b*x, X, Y, x);

where X is your list or array of x data and Y is your list or array of y data. If you're using CurveFitting:-LeastSquares, then the command is

CurveFitting:-LeastSquares(X, Y, x, curve= b*x);

Surely this problem was contrived to point out the flaws in a floating-point arithmetic system. In particular, this problem looks contrived specifically to dupe Maple because the error occurs right at 10 decimal places, which is Maple's default. So I can't believe that it is truly a BIG problem for you.

The problem occurs just in the computation of

.5 - 1234567892;

This value requires 11 digits to represent, obviously, but you only have 10 digits to work with (by default). It is easy to increase the default:

Digits:= 11;

Now you will get accurate results for this computation.

Another option is to use exact arithmetic:

convert(.1234567891, rational, exact)*10^10 +1/2 - 1234567892;

Exact arithmetic works accurately regardless of the value of Digits.

I managed to make it into a simulation for any number of doors and still simplify and speed up the code some.

# Proc to make a uniform random selection from a set.
Rand:= (S::set)-> `if`(nops(S) < 2, S[], S[(rand() mod nops(S)) + 1]):

# Simulation for the n-door Monty Hall problem.
# Returns the number of trials the player wins.
MontyHall:= proc(
     {trials::posint:= 10000},
     {doors::posint:= 3},
     {switches::list(truefalse):= [true]},
     $
)
local
     k, wins:= 0, Doors,
     PrizeDoor, #The door with the prize behind it.
     PlayersDoor,  # The door the player picks.
     AllDoors:= {$1..doors}   
;
     to trials do
          Doors:= AllDoors;
          PrizeDoor:= Rand(Doors);
          PlayersDoor:= Rand(Doors);
          for k to doors-2 do
               # Monty opens a door, which is thus removed from the set
               # of available doors.
               Doors:= Doors minus {Rand(Doors minus {PrizeDoor,PlayersDoor})};       
               if switches[k] then
                    PlayersDoor:= Rand(Doors minus {PlayersDoor})
               end if
          end do;
          if PlayersDoor = PrizeDoor then  wins:= wins+1  end if
     end do;
     wins
end proc:     

MontyHall(trials= 100000, doors= 5, switches= [false,false,true]);
                             79904

For longer sections of code within a worksheet, use Code edit Regions. These are available from the Insert menu. They do allow the tab key. They offer syntax-based color coding and highighting and bracket matching.

One drawback: In M17, Code Edit Regions do not place the cursor near a syntax error when you make one; you're on your own for finding syntax errors. This is a major drawback that M16 did not have.

To excecute a Code Edit Region, use Ctrl+E. The option is also accessible from the context menu.

Here is my solution to the 4-door Monty Hall problem. I managed to simplify the code a lot compared to my 3-door solution.

# Proc to make a random selection from a set.
Rand:= (S::set)-> S[(rand() mod nops(S)) + 1]:

# Simulation for the 4-door Monty Hall problem. N is the number
# of trials. Returns the number of trials the player wins.
MontyHall4:= proc(
     N::posint,
     {switch1st::truefalse:= false},
     {switch2nd::truefalse:= false},
     $
)
local
     k, wins:= 0, Doors,
     PrizeDoor, #The door with the prize behind it.
     PlayersDoor,  # The door the player picks.
     AllDoors:= {$1..4}   
;
     to N do
          Doors:= AllDoors;
          PrizeDoor:= Rand(Doors);
          PlayersDoor:= Rand(Doors);
          # Monty opens a door, which is thus removed from the set
          # of available doors.
          Doors:= Doors minus {Rand(Doors minus {PrizeDoor,PlayersDoor})};       
          if switch1st then
               PlayersDoor:= Rand(Doors minus {PlayersDoor})
          end if;
          # Monty opens another door:
          Doors:= Doors minus {Rand(Doors minus {PrizeDoor,PlayersDoor})};
          if switch2nd then
               PlayersDoor:= Rand(Doors minus {PlayersDoor})
          end if;
          if PlayersDoor = PrizeDoor then  wins:= wins+1  end if
     end do;
     wins
end proc:     

MontyHall4(100000);
                             25092
MontyHall4(100000, switch1st);
                             37562
MontyHall4(100000, switch1st, switch2nd);
                             62565
MontyHall4(100000, switch2nd);
                             74751

If the blocks of a partition cannot be made to have equal products, then we seek to minimize the "variation" (defined below). I will use Brian's example set

S:= {3, 4, 5, 6, 8, 9, 28, 30, 35}.

First, I compute the ideal target product for a block---the cube root of the product of the elements of S. Call this p. Then for each block b, I measure its deviation from this ideal as abs(ln(`*`(b[]) - ln(p)). I use ln because I think that being a factor of x too large should be considered the same deviation as being a factor of x too small. For each partition P, I measure the variation as the sum of the deviations of its blocks.

The partitions are generated with Joe's Iterator package. You'll need to download this from the Maple Applications Center if you don't already have it.

restart:
S:= {3, 4, 5, 6, 8, 9, 28, 30, 35}:
AllP:= [seq(P, P= Iterator:-SetPartitions(S, [[3,3]], compile= false))]:
lnp:= evalf(ln(`*`(S[])^(1/3))):
Var:= proc(P::list(set))
local r:= evalf(`+`(map(b-> abs(ln(`*`(b[]))-lnp), P)[]));
end proc:
sort(AllP, (x,y)-> Var(x) < Var(y))[1];
              [{3, 9, 35}, {4, 8, 28}, {5, 6, 30}]
---which agrees with Brian's result.

Rather than using sort, I wanted to use Joe's trick of setting the attribute of a float and then using min. But this does not work:

Var:= proc(P::list(set))
local r:= evalf(`+`(map(b-> abs(ln(`*`(b[]))-lnp), P)[]));
     setattribute(r,P)
end proc:

The procedure returns garbage which displays as a question mark. I am hoping that Joe can tell me why.

@fluff

Forget about the identical(0,1,2) thing. It won't work, because I didn't realize that you allowed the possibility of the player not switching in the first round and then switching in the second. You need two switch parameters:

MontyHall:= proc(N::posint, {switch1st::truefalse:= false}, {switch2nd::truefalse:= false})

Then is the code you'll have

if switch1st then
. . .

if switch2nd then
. . .

Also, you still need to correct rendomize to randomize, as pointed out by Acer.

It looks like the command defining params is not being executed. Try placing the restart in its own execution group.

Your newsys is a list. In order to use it with union, it must be a set. In your dsolve command change newsys to {newsys[]}.

Download Leontief2.mw

Download odesys.mw

GraphTheory:-Graph will accept an adjacency matrix but not an incidence matrix. Here is a procedure that will convert an incidence matrix into an adjacency matrix.

IncidenceToAdjacency:= proc(J::Matrix)
uses LA= LinearAlgebra;
local
     r:= LA:-RowDimension(J),
     A:= Matrix(r,r, shape= symmetric),
     c:= LA:-ColumnDimension(J),
     j, ones
;
     for j to c do
          ones:= ArrayTools:-SearchArray(J[..,j]);
          A[ones[1], ones[2]]:= 1
     end do;
     A
end proc:

(This procedure only handles the case of an unweighted, undirected graph.)

Using your example incidence Matrix:

M:=Matrix([[1,1,0,1,0,0,0],[0,1,1,0,1,0,0],[1,0,1,0,0,1,1],[0,0,0,1,1,1,0],[0,0,0,0,0,0,1]]):
A:= IncidenceToAdjacency(M);

G:= GraphTheory:-Graph(A):
GraphTheory:-DrawGraph(G);

Here is a complete simulation for the three-door Monty Hall problem (the classic). The code is very straightforward to read. See if you can extend this to the four-door problem.

MontyHall:= proc(N::posint, {switch::truefalse:= false})
# Simulation for the 3-door Monty Hall problem. N is the number
# of trials. Returns the number of trials the player wins.
local
     wins:= 0,
     PrizeDoor, #The door with the prize behind it (1..3).
     MontysDoor,  # The door that Monty reveals (1..3).
     Players1stDoor,  # The door the player picks first.
     Players2ndDoor,  # The door the player picks second.
     Doors:= {1,2,3},
     rand3:= rand(1..3), rand2:= rand(1..2)
;
     to N do
          PrizeDoor:= rand3();
          Players1stDoor:= rand3();
          if Players1stDoor = PrizeDoor then
               MontysDoor:= (Doors minus {PrizeDoor})[rand2()]
          else
               MontysDoor:= (Doors minus {PrizeDoor,Players1stDoor})[]
          end if;
          if switch then
               Players2ndDoor:= (Doors minus {MontysDoor,Players1stDoor})[]
          else
               Players2ndDoor:= Players1stDoor
          end if;
          if Players2ndDoor = PrizeDoor then
               wins:= wins+1
          end if
     end do;
     wins
end proc:     

MontyHall(100000);
                             33282
MontyHall(100000, switch);
                             66915

You asked:

How would i label the axis and change its font and size?

By using plot options labels and labelfont. (Example shown below.)

Also, how would i move the horizontal axis down, so that the whole diagram can be seen?

By using plot options axes= box and view. (Example shown below.)

Finally, there are some odd points between 0 and 1 on the vertical axis. How do you get rid of these?

Those are very interesting and seem to be an integral result of the computation for small values of r and some x, although it may also be an effect of round-off error. I will investigate that further. To eliminate them from the plot, slightly increase the lower bound for the view of the horizontal axis. I used 0.1 as the increase amount in the example show below.

Bonus: Here's how you can generate this plot in well under 1 second by using evalhf. It's at least a factor-of-50 improvement timewise over the old way. I set this up such that it will be easy to modify for other bifurcation diagrams.

restart:
Rloop:= proc(Pts::Matrix, r_min::realcons, r_max::realcons, N::posint, M::posint)
local n, r, x, r_range:= r_max - r_min;
     for n to N do
          r:= r_min+n/N*r_range;
          x:= Pts[n,2];
          to M do  x:= x*exp(r*(1-x))  end do;
          Pts[n,1]:= r;  Pts[n,2]:= x
     end do;
     NULL
end proc:     

(N,M):= 10^~(4,2):  (x_min,x_max):= (0,1):  (r_min,r_max):= (0,4):
bifpoint:= Matrix(N,2, datatype= float[8]):
bifpoint[.., 2]:=
     < RandomTools:-Generate(list(float(range= x_min..x_max, method= uniform), N))[] >:
evalhf(Rloop(bifpoint, r_min, r_max, N, M)):
pitchf:= plots:-pointplot(bifpoint, symbol=point):
plots:-display(
     pitchf,
     axes= box, view= [r_min+0.1..r_max, -1..5],
     labels= ['r','x'], labelfont= [TIMES,BOLDITALIC,16],
     axesfont= [HELVETICA,BOLD,12]
);

You can use the events option to dsolve to detect when the derivative is 0. If you post the equations, I'll help to set it up. Otherwise, you can read about it at ?dsolve,events . (Warning: This help page is quite esoteric.)

While it's true that you can't use solve with the output of dsolve(..., numeric), you can use fsolve.