Hot answers tagged

29

Usage Just use this function with any polyhedron in in form: GraphicsComplex[pts_, Polygon[vertices_, ___]]. When I find time and motivation maybe I will add more DownValues so it can be more general. At the moment you can play with solids given by PolyhedronData[... "Faces"]: polyhedronRandomWalk[ PolyhedronData["DuerersSolid", "Faces"] ] ...


23

BernoulliDistribution is a perfect fit for this. RandomVariate[BernoulliDistribution[1 - 0.1], {50}] {1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1} Also, as kguler states, you can use RandomChoice, but the benefit of BernoulliDistribution is that ...


12

ClearAll[dist] dist [rho_] := CopulaDistribution[{"Binormal", rho}, {UniformDistribution[{0, 1}], UniformDistribution[{0, 1}]}]; data1 = RandomVariate[dist[-.9], 5000]; ListPlot[data1, Frame -> True, AspectRatio -> 1] data2 = RandomVariate[dist[.6], 5000]; ListPlot[data2, Frame -> True, AspectRatio -> 1]


12

Your code shows you are somewhat confused about the indices. Here is some code that is more Mathematica idiomatic and which makes keeping the indices straight much easier. I am running the simulation with a 24 hour step to cut down the data plotted. Should work just a well with your time step. maxk = 10; mink = 2; steps = 365; batterylevel = ...


11

In Mathematica it is natural to approach such a task with list operations and pattern matching. dice1 = RandomInteger[{1, 6}, 100]; dice2 = RandomInteger[{1, 6}, 100]; Count[dice1 + dice2, 2 | 3 | 4 | 5 | 6] You seem to be a very new beginner, since you are using x[i] and y[i] as if these are vectors, when they are in fact not, in Mathematica. Mathematica ...


11

There are many ways you can do this, e.g. ri = RandomInteger[{1, 6}, {100, 2}]; SortBy[Normal@GroupBy[ri, Total, Length@#/100. &], First] yielding: {2 -> 0.01, 3 -> 0.05, 4 -> 0.11, 5 -> 0.08, 6 -> 0.13, 7 -> 0.12, 8 -> 0.17, 9 -> 0.14, 10 -> 0.11, 11 -> 0.07, 12 -> 0.01} rules linking sum to frequency. You can also exploit ...


7

The Finite Element solver in Mathematica does run in parallel, both element computation and the linear solve process are spread over the CPU cores available. Additionally, the option "MeshElementBlocks" for ToElementMesh splits the mesh elements in blocks which could be used for a domain decomposition. To get a more detailed answer you'd need to clarify ...


7

While the answers so far have covered a lot of ground already I have not seen EmpiricalDistribution. I would like to build upon this observation by providing a couple of general considerations that I have found to be useful when doing statistical experiments using Mathematica. What users of Mathematica may take for granted may surprise newcomers: You can ...


7

When building a simulation like yours you should test the performance of the individual components before incorporating them into the simulation. That is, you should know the cost of the components as well as the values they return. Here is a simple example based on your code. You use Clip in a couple of places to limit values from below. That is suspect, ...


7

Here is a modified version of your code. On my PC it completes your example run in about 5 seconds. I won't try to describe every change but will point out the major features. Some of the changes are stylistic rather than performance-based. This is not a criticism of your style but a reflection of the way I broke the original code down in order to ...


6

You could use Region functionality (for simpler regions), e.g. rw[pt_, s_, n_, reg_] := Module[{ch = {{0, 0}, {1, 0}, {-1, 0}, {0, 1}, {0, -1}}, np, st}, st = RandomChoice[ch, n]; FoldList[If[RegionMember[reg, #1 + s #2], #1 + s #2, #1] &, pt, st] ] an[p_, step_, num_, regn_] := With[{pnts = rw[p, step, num, regn]}, ListAnimate[ ...


5

Let me give a model solution which can easily be adapted. 2 dimensions Consider a photon which moves in the x-y-plane, starting at time t = 0 in the origin {0,0} and moving towards the positive x-axis. At each tick of the clock, corresponding to a constant distance 1 travelled by the photon, the photon will experience a scattering event which leads to a ...


5

My interpretation is: the photon will be scattered by an angle $\alpha$ (given by getScatterAngle), and the deviation will occur with equal probability in every direction. (For example, a photon initially going along the $z$ axis will be rotated by an angle $\alpha$ about a randomly chosen axis that lies in the $x,y$ plane.) When I've written Monte Carlo ...


5

Take a uniform random distribution and check if it is above some threshold (0.9 in your case). For example: dist[] := If[RandomReal[] > 0.9, 0, 1]; Table[dist[], {i, 100}]


5

First, you're not using a fixed step method. (An Euler scheme may be applied to any step size and to one that varies.) To get a true fixed step method you have to turn off "DiscontinuityProcessing" when you have a discontinuous ODE; otherwise, NDSolve will try to adapt the steps to account for the discontinuity. The "DiscontinuityProcessing" stage resets ...


5

You go to the casino. You bet a dollar the ball will land on red. If it does you get two dollars if it does not you get zero. Likewise if you bet the ball will land on black. You cannot bet it will land on green. Use this previous answer to help do the pattern matching. Note to the original poster, I found that using Google and searching for mathematica ...


4

I misunderstood and commented only about computing the Cesaro means as per the question In s I have all partial sums, but I do not know how to divide them by corresponding n. The desired scatter plot of 500 Monte Carlo attempts with samples of increasing length could be obtained with something like ticks = Range[-0.06, 0.06, 0.02]; s = Table[u = ...


4

In your example plot, the Monte Carlo integral is computed afresh for each new amount of sampling points: s[n_] := Total[Sqrt[1 - #^2] & /@ RandomReal[1, n]]/n Then take one random point, find the mean, take two new random points, find their mean, and so on until 500. The result is: ListPlot[Table[s[n], {n, 500}], PlotRange -> {0.7, 0.9}] The ...


4

You can use RandomVariate to sample from a DiscreteUniformDistribution and then add up the pairs, calculate the probability of the sums observed, and then extract the probabilities of interest. (#/100. & /@ Counts[Plus @@@ RandomVariate[ DiscreteUniformDistribution[{1, 6}], {100, 2}]] )[#] & /@ Range[2, 6] Hope this helps.


4

The reason it's not finishing is that you set MaxSteps -> Infinity with a high AccuracyGoal. We can actually solve this system analytically, by replacing NDSolve with DSolve: qq = DSolve[{q'[t] == v[t]/r - q[t]/r*(1/c), q1'[t] == (((q[t] + q1[t])/(2*a*ϵ0*k2))*δ), q[0] == 10*10^-9, q1[0] == 10*10^-9}, {q, q1}, t]; We can inspect the solution: ...


4

RandomChoice[{.9, .1} -> {1, 0}, 10] (* {0, 1, 1, 1, 1, 1, 1, 0, 1, 1} *) Timing results Timing[RandomVariate[BernoulliDistribution[.9], {10^8}];] (* {3.38014, Null} *) Timing[RandomChoice[{.9, .1} -> {1, 0}, 10^8];](* {5.64937, Null} *) dist[] := If[RandomReal[] > 0.9, 0, 1]; Timing[Table[dist[], {i, 10^8}];] (* {13.9842, Null} *) One nice ...


4

dist = TransformedDistribution[x + y, { Distributed[x, DiscreteUniformDistribution[{1, 6}]], Distributed[y, DiscreteUniformDistribution[{1, 6}]]}]; SeedRandom[1] For small sample sizes, the match to the theoretical values is poor. data = Total /@ RandomInteger[{1, 6}, {100, 2}]; Show[ Histogram[data, {1.5, 12.5, 1}, "PDF", ...


3

data10 = RandomVariate[f[13, 0.5], {10, 25}]; (* 10 data sets from distribution f *) lls= LogLikelihood[EstimatedDistribution[#, f[a, b], {{a, 1}, {b, 1}}], #] & /@data10 (*{-32.4994, -25.2268, -21.9671, -26.8963, -25.9164, -22.8958, -26.5247, -24.9622, -33.9319, -28.6512}*) maxll=Block[{k=1}, MaximalBy[Last][ With[{dist = ...


3

I'll try it as an answer though it's short... getTrajectory[startV_, steps_] := Accumulate@NestWhileList[stepVector, startV, RandomReal[] >= 1/steps &];


3

You can also use RandomVariate with DiscreteUniformDistribution: rW[a_, b_, n_] := Accumulate[Prepend[ RandomVariate[DiscreteUniformDistribution[{{-a, a}, {-b, b}}], n]], {0,0}] dt = rW[10, 20, 100]; Graphics[{PointSize[Large], Red, Point@#, Thick, Blue, Line@#} &@dt, Frame -> True, Axes->True, AspectRatio -> 1/GoldenRatio] We get the ...


3

f[r_?NumberQ,n_Integer]:={First[#],#.{r,Sqrt[1-r^2]}}&/@RandomReal[NormalDistribution[0,1],{n,2}]; Produces n pairs of numbers with the correlation r.


3

nextGen[n_] := Total@RandomChoice[{11/32, 3/8, 3/16, 3/32} -> {0, 1, 2, 3}, n] simulate[n0_, nrOfGenerations_] := Total@NestList[nextGen, n0, nrOfGenerations] Now we can simulate six generations a hundred times and compute the mean value. The initial number of organisms is 10 in this example. Table[simulate[10, 6], {100}] // Mean // N (* Out: 75.42 *) ...


3

Here's a q-d implementation: With[{hits = Cases[Partition[RandomChoice[{17/36, 17/36, 1/18} -> {1, 2, 3}, 1000000], 5, 1], {x_, x_, x_, x_, y_}]}, Count[hits, {x_, x_, x_, x_, y_} /; y != x] - Length@hits] It will run through 1 million spins, checking for any 4-runs and counting those where the subsequent spin differs, less the count in total (so ...


2

What you are missing, I believe, is sufficient experience of Mathematica's core language at the functional level that you experimenting with. I give you credit for making a good try at formulating your code in a functional way, but I'm afraid you gone somewhat wide of the mark. I have put your code into a form that works and which I think preserves your ...


2

Maybe you can try a matrix approach. 1/ The idea is to generate a matrix like this one : mat = {{a, 1, 0}, {b, 0, 1}, {c, 0, 0}}; then you can see that: {x1, x2, x3}.{{a, 1, 0}, {b, 0, 1}, {c, 0, 0}} {a x1 + b x2 + c x3, x1, x2} gives you the format you want. 2/ Then you can use that matrix directly in NestList (without the need to define a ...



Only top voted, non community-wiki answers of a minimum length are eligible