数値的 Mathematica
Numerical Mathematica

メインテーマ:
Accuracy（確度）とPrecision（精度）の使用法
Main Theme:
Using Accuracy and Precision

Mathematica の数値コマンドを用いると結果のクォリティを詳細に設定することがでますが，通常はAutomatic設定で十分です．Automatic設定が理想的でない状況の1つに，2つ以上の数値コマンドを同時に用いたときが挙げられます．例えば数値コマンドに関数値を提供するために他の数値コマンドを用いるときなどです．これらのアプリケーションを効率よく使用するために，Mathematica におけるaccuracy及びprecisionの意味と，これらがAccuracyGoalやPrecisionGoalオプションとどのように関わってくるかを理解する必要があります．ここではまずMathematica のaccuracy及びprecisionについて紹介し，緻密な誤差制御を行うためにこの概念が用いられる例を2つ示します．

Mathematica's numerical commands allow excellent control over result quality and usually the Automatic settings are sufficient for most uses.  One situation which arises where the Automatic settings are not ideal is when you are using more than one numerical command at a time -- e.g. using one to provide a function value for another.   For these types of applications to be efficient, you need to have a good understanding of the meaning of now Mathematica interprests accuracy and precision and how these are related to the AccuracyGoal and PrecisionGoal options.  In this talk, I will guide you first through an introduction to accuracy and precision in Mathematica and then two examples where the concepts are used for tight error control.

初期化
Initialization

Off[General::spell];
Off[General::spell1];

$MinPrecision  =  -Infinity;

$Path  =  
        Join[$Path,  {"C:\\Documents and Settings\\rknapp\\My \
Documents\\Talks\\2001 IMS", "C:\\Documents and Settings\\rknapp\\My \
Documents\\Talks\\2001 IMS\\Packages"}];

Needs["Utilities`Scale`"]

accuracyとprecisionの基礎
Basics of Accuracy and Precision

まず始めにMathematica ブックの関連セクション3.1.4をご覧ください．

A starting point is to look at the relevant section from The Mathematica Book: Section 3.1.4.

これらの2つの関係を理解するためにはscaleの概念が便利です．数値のscaleは小数点の左側にある数字の数です．そのため，実軸上のあらゆる数値に関して以下の式が成り立ちます．

The concept becomes clearer if you use the concept of scale to form a relationship between the two.
The scale of a number is the number of digits it has to the left of the decimal point.  Thus, for any number on the real line x, we have

Precision[x] == Scale[x] + Accuracy[x]

以下に数量を図解します．

Here is a pictorial representation of the quantities:

(Overscript[(Overscript[(Underscript[\(\(\(\ \)\(x\_1\ x\_2\ ⋯\ x\_s\)\)\+⏟\), Scale]).(Underscript[\(\(x\_\(s + 1\)\ x\_\(s + 2\)\ ⋯\ x\_\(s + a\)\)\+⏟\), Accuracy]), ⏞]), Precision])

Accuracy，Precision，及びScaleのいずれも負の値をとることができます．

Any of the values for Accuracy, Precision and Scale can be negative.

これらの関係を，いくつかの例を通して見てみましょう．

Lets see how this relationship holds for some examples

PSA[x_]  :=  {Precision[x],  Scale[x],  Accuracy[x],  
        Precision[x]  == (Scale[x]  +   Accuracy[x])}

このMathematica ブックの例

The Mathematica Book example ...

PSA[xacc = 431.123145333555141444]

はほぼ正しく動作するのですが，問題は，PrecisionとAccuracyは整数値ではなく，おおよその数量として考えるべきだという点です．

... almost works.  The problem is that you really need to consider Precision and Accuracy as approximate quantities, not integral numbers of digits:

PSA[x_]  :=  {Precision[x,  Round -> False],  Scale[x],  
        Accuracy[x,  Round -> False],  
        Precision[x,  
            Round -> 
              False]  == (Scale[x]  +   Accuracy[x,  
                Round -> False])}

今度は

and now...

PSA[xacc]

正しい結果を表示します．

... comes out as it should!

これらは厳密な値でも計算可能です．

This even works for exact numbers

PSA[1]

小さい値 (<< 1)では，scaleは負であることに注目してください．

Notice that for small (<< 1) numbers, the scale is negative

PSA[10. ^-1000]

これはprecisionとaccuracyがある特定の桁数であると考えるより，浮動小数点数であると考えた方が妥当であるからです．ある特定のaccuracyをもつ数は，その数を中心とするインターバルを表していると考えることができます．

Since it clearly makes more sense to think of Precision and Accuracy as floating point quantities than it does as a particular number of digits.  A  number with a particular accuracy can best be thought of as representing an interval which is centered at the number.

数値xの属性がPrecision[x] == p及びAccuracy[x] == aである場合を考えてみます．するとxの誤差の絶対幅は(10^(-a))になります．つまり，xは以下の誤差範囲中のあらゆる\(TraditionalForm\`\(x\_*\)\)を表します．

Suppose that x has Precision[x] == p  and Accuracy[x] == a.  Then the absolute width of the error interval for x is (10^(-a)).  This means that x represents any number \(TraditionalForm\`\(x\_*\)\)in the interval

\(TraditionalForm\`x - \(1\/2\) 
          10\^\(-a\) ≤ \(x\_*\) ≤ 
      x + \(1\/2\) 10\^\(-a\)\)

または，同様に以下の式も成り立ちます．

or equivalently

\(TraditionalForm\`x(
        1 - \(1\/2\) 10\^\(-p\)) ≤ \(x\_*\) ≤ 
      x(1 + \(1\/2\) 10\^\(-p\))\)

これを具体的にするために，上で定義されたxaccを考えます．ここで誤差の幅は以下の通りであるとします．

To make this concrete, consider the number xacc defined above, where the error interval has width

10^-Accuracy[xacc,  Round -> False]

ですから，xaccは以下の範囲内のどの\(TraditionalForm\`\(x\_*\)\)にもなり得ます．

and so xacc could be any number \(TraditionalForm\`\(x\_*\)\) in

xinf  =  
        SetPrecision[xacc,  Infinity];
iinf  =  
        SetPrecision[10^-Accuracy[xacc],  Infinity]/
          2;
\(xinf\  - iinf\  ≤ \ \(x\_*\) ≤ \ 
      xinf\  + \ iinf\)

誤差の境界をシャープで明確にするために有理数を用います．

Rational numbers have been used here to make the boundary of the interval sharp and clear.

ここで負のprecisionの概念を例証します．まず，以下のものを考えてみてください．

Here is an example which illuminates the concept of negative precision.  Consider

poorx  =  3`-5

これは以下の誤差範囲中のあらゆる数になり得ます．

This could be any number in the interval

pinf  =  
        SetPrecision[SetPrecision[poorx,  100],  
          Infinity];
pint  = 10^-Precision[poorx];
\(pinf*\((1\  - pint/2)\)\  ≤ \ \(x\_*\) ≤ \ 
      pinf*\((1\  + \ pint/2)\)\)

誤差の大きさはaccuracyと一致します．

The error interval size is consistent with the accuracy

150003  +  149997  ==  
      10^-Accuracy[poorx,  Round -> False]

precisionの概念はゼロ付近で無効化するものの，accuracyは絶対で，どこでも使うことができます．

While the idea of precision breaks down around zero, accuracy is absolute so it can be used anywhere.

これらの考え方や，Mathematica において有効算法を形成するためにこれらがどのように使われているかをもっとお知りになりたい方は，1998 Chicago conferenceにおけるMark Sofroniuの講演「Efficient and Numerically Stable Computations」をご参照ください．

If you would like to learn more about these ideas and how they are used to form our significance arithmetic, and excellent resource is Mark Sofroniou's talk on Efficient and Numerically Stable Computations from the 1998 Chicago conference.

PrecisionGoalとAccuracyGoal
PrecisionGoal and AccuracyGoal

accuracyとprecisionの概念は，AccuracyGoal及びPrecisionGoalオプションを用いて数値関数まで拡張できます．以下にこれを解明する例を示します．

The concepts of accuracy and precision are extended to numerical functions with options AccuracyGoal and PrecisionGoal.  An example helps illuminate this

NIntegrate[Sin[2 π x], {x, 0, 1}]

デフォルトの基準を満たしていないため，警告を表示します．

Gives a warning message because it could not meet its default criteria:

Options[NIntegrate]

この結果はとてもよいように思われるかもしれませんが，相対誤差に関してはよくありません．ここに適した基準は絶対誤差なので，以下のものを用いることができます．

The result may seem quite good, but in terms of relative error it is poor.  The appropriate criterion to use here is absolute error, so you could use

NIntegrate[Sin[2 π x], {x, 0, 1},  
      AccuracyGoal -> 6]

一般的に，PrecisionGoalとAccuracyGoalの両方が使われている場合，2つのうち厳しくない方の基準が満たされれば，その結果は許容範囲とみなされます．

And we are fine.  In general, when both Precision and AccuracyGoal are present, the result is deemed OK if the least stringent of the two is estimated to be met.

通常，Precisionは以下の

In general, Precision is the better criterion, since it allows both

NIntegrate[Exp[x], {x, 100,  101} ]

及び

and

NIntegrate[Exp[x], {x, -101,  -100}]

の両方が可能なので，同様の相対（もちろん絶対ではありませんが）誤差に対して妥当な解を得るためによりよい基準となります．

To give reasonable answers with the same relative (though certainly not absolute!) error.

簡単な例: 割り算
Simple Example: Division

単純な割り算は，掛け算と足し算を用いてニュートン法のみで行われます．

Simple division can be done using multiplication and addition only by using Newton's method.  With

f[a_][x_]  :=  a  -  (1/x)

であるすると，

then

Expand[x  -  f[a][x]/f[a]'[x]]

であるので，以下のように定義できます．

so you can define

DivisionStep[a_][x_]  := (2  - a x) x

ここで，固定小数点を見つけることができます．

Then you can look for the fixed point

pin  =  N[Pi];
FixedPoint[DivisionStep[pin],   .1]

上式は機械精度の数値で始まっているので，ここの計算には機械精度が用いられます．bignumで始まる場合は，より精度の高い結果が得られます．

Machine numbers were used in the computation above because it was started with a machine number.  If you start with a bignum, then you can get a higher precision result.

prec  =  100000;
pin  =  N[Pi,  prec];
start  =  N[1/10, prec];
Timing[rip  =  
          FixedPoint[DivisionStep[pin],  start];]

これは少々時間がかかりますが，以下の精度の高い結果を返します．

which takes a little while, but does give a result which is close

rip*pin  -  1

ニュートン法は，反復演算が収束すると誤差が消滅するという意味で自己修正アルゴリズムです．したがって，低精度から始めて，根に近付くに従って精度を高めることが可能です．

Newton's method is a self correcting algorithm in the sense that if the iterations are converging, errors will disappear.  Thus, it is possible to start with lower precision and increase the precision as you get closer to the root.

AdaptiveNewton[step_,  start_,  
        opts___]  :=  
Module[{goal,  
          prec}, 
c  =  0; 
        goal  =  
          PrecisionGoal  /.  
            Flatten[{opts,  
                Options[AdaptiveNewton]}]; 
prec  =  
          Min[Precision[start],  goal]; 
SetPrecision[
          NestWhile[AdaptiveNewtonStep[step,  prec],  start,  
            Function[prec  <  2*goal]],  
          goal]
]

FixedPointを使うこともできますが，理想的ではありません．なぜならその停止のための判定条件は，結果の有効数字が変化していないということだからです．ここでの考え方は，精度を高めて有効桁数を増やすというものなので，プロセスは停止しません．NestWhileは停止条件をテストするために関数に適用できるため，非常に適しています．これを用いて，各々の反復ステップを行う関数を定義する必要があります．

FixedPoint could be used, but it is not ideal because its criterion for stopping is that the significant digits of the result is not changing. Since the idea is to increase the precision and add significant digits, the process would never terminate. NestWhile fits the bill very nicely because it allows you to put in a function to test for a stopping condition.  Using it, we need to define a function which handles each iterative step.

\(\(AdaptiveNewtonStep[step_, \ prec_]\)[\ x_]\  := \ 
      Module[{xnew, \ Δx}, 
xnew\  = \ 
          step[x]; 
Δx\  = \ 
          xnew\  - \ 
            x; 
If[\((Δx\  == \ 
                0)\) || \ \((prec\  < \ 
                2 \((Scale[
                        x] - \ \ Scale[Δx])\))\), 
\ 
          prec\  *= \ 
            2. ; 
print["Increasing Precision to ", \ 
            prec, \ " at iteration ", \ \(++c\)]; 
xnew\  = \ SetPrecision[xnew, \ prec]]; 
xnew]\)

ここでは2つの便利な方法を用います．

There two useful tricks being used here.

引数の参照 
Arguments by reference

まず，precをTimesByで変化させます．precは記号でなければなりません（また，値をもたなければなりません．そうでないと再帰誤差が生じます）．以下の計算を行う場合，

First, we are changing prec with TimesBy.  This requires that prec be a symbol (and also that it has a value, since if not a recursion error will occur).  But if you try it

prec  =  $MachinePrecision; {AdaptiveNewtonStep[
          DivisionStep[pin],  prec][1. /Pi], prec}

デフォルトとして引数が評価されるため，上記のような結果になります．Unevaluatedを用いることもできると思われるかもしれませんが，以下の式

That is because by default arguments are evaluated.  You might expect that it would be possible to use Unevaluated, but this gives

{AdaptiveNewtonStep[DivisionStep[pin],  Unevaluated[prec]][
        1. /Pi], prec}

ではUnevaluatedを用いるのは難しくなります．関数をどのように用いた場合でもprecは評価されるべきではないので，これを関数自体に関連付けた方がよいと考えられます．

Using Unevaluated can be difficult at best.  Since any use of the function requires that prec not be evaluated, it is better to associate this with the function itself anyway.

SetAttributes[AdaptiveNewtonStep,  HoldAll];

これで期待通りの計算が行われます．

Now everything works just as we want.

{AdaptiveNewtonStep[DivisionStep[pin],  prec][1. /Pi], 
      prec}

このステップがまず最初に評価されるべきである複雑な関数であったとしても，それはAdaptiveNewtonを用いたときにすでに評価されるので，問題ありません．

Even if the step were a complicated function which should be evaluated first, that is not a problem since it will have already been evaluated when using AdaptiveNewton

Mathematica ではデフォルトでは引数は数値で受け渡されます．例えば，

In Mathematica, the default is to pass arguments by value.  For example

Subtract[Times[rip,  Pi],  1]

は再帰的に式を減少させ，引数を評価し，それから関数を評価することにより計算されます．Times及びSubtractのコードはrip，Pi，及び1が変数と関連付けられているかどうかを知りません（あるいは必要ありません）．

is computed by descending the expression recursively, evaluating arguments, then evaluating the functions.  The code for Times and Subtract have no (nor do they need to) knowledge about whether rip, Pi, and 1 are associated with variables.

また，一方で

On the other hand,

Set[a,  5]

ではaが記号であることが必要です．aにはすでに値があるかもしれませんが，ここで5が代入されます．Setは記号aを設定します．すると設定後，aは数値5の参照となります．Mathematica では参照は記号と関連付けられています．

requires that a be a symbol.  a could have a previous value, which will be replaced by 5.  Set works with the symbol a, which after set is run is a reference to the number 5.  In Mathematica references are associated with symbols.

また，

What about

b[1]  =  a

では，記号ではないb[1]に割り当てているように見えます．実際にはルールはbに付属しています．

In this case it appears that we are assigning to b[1], which is not a symbol.  In reality, the rule is attached to b:

DownValues[b]

しかし，Mathematica では評価を遅らせることによって参照を渡すことができ，とても便利であるということにご注目ください．

However, it is important to note that it is possible to pass by reference in Mathematica by delaying evaluation and this is a useful technique to be aware of.

おおよそのscaleを用いる
Using approximate scale

不等式

The inequality,

prec  <  2 (Scale[x] -   Scale[Δx])

は極めて曖昧ですが，scaleの定義を用い，

probably seems quite obscure, but if you use the definition of scale,

prec  <  
      2 (Log[10. ,  Abs[N[x]]]  -   Log[10. ,  
              Abs[N[Δx]]])

さらに式の両辺を-2で割って指数をとると以下のようになります．

and you divide both sides of the equation by -2 and exponentiate, this becomes

Abs[Δx]/( Abs[x]) <  
      (10^(-prec/2))

左辺はxに対するΔxの大きさであり，右辺はprecに使用可能な桁数の半分の（accuracyの）許容誤差範囲を表しています．

The left-hand side is the size of Δx relative to x and the right hand side represents an (accuracy) tolerance of half of the digits available in prec.

ここで使用可能な桁数を半分にしたことは，精度を倍にすることに関連しており，ニュートン法の二次方程式の収束に直接関係があります．(Subscript[x, 1]), (Subscript[x, 2]), ...が根の極限値\(TraditionalForm\`\(x\_*\)\)に収束するニュートン反復列であるとすると，以下のようになります．

Why half?  This is tied in with our doubling of the precision, and is associated directly with the quadratic convergence of Newton's method: If (Subscript[x, 1]), (Subscript[x, 2]), ...  are a sequence of Newton iterates converging to a root limit \(TraditionalForm\`\(x\_*\)\), then

\[LeftBracketingBar](Subscript[x, k + 1]) - 
          (Subscript[x, k])\[RightBracketingBar] -> \[LeftBracketingBar](Subscript[x, k]) - (Subscript[x, k - 1])(\[RightBracketingBar]^2)

したがって，Δ(Subscript[x, k]) = |(Subscript[x, k]) - 
          (Subscript[x, k - 1])|が(10^(-prec/2))（つまり差によってキャンセルされる最初のprec/2）の位数であるなら，次のΔ(Subscript[x, k + 1])= |(Subscript[x, k + 1]) - 
          (Subscript[x, k])|は(10^(-prec))の位数であると想定されます．つまり(Subscript[x, k + 1])が(Subscript[x, k])以外の有効数値をもつとすると，精度を高めた方がよいということになります．

Thus, if Δ(Subscript[x, k]) = |(Subscript[x, k]) - 
          (Subscript[x, k - 1])| is of order (10^(-prec/2))(i.e. the first prec/2 digits cancel in the difference) you expect the next Δ(Subscript[x, k + 1])= |(Subscript[x, k + 1]) - 
          (Subscript[x, k])| to be of the order of (10^(-prec)), so if you want (Subscript[x, k + 1])to have any significant digits different than (Subscript[x, k]), you had better increase the precision!

現在の全精度で余分な計算を行うことを回避するためにscaleを用います．例えば上で導いた不等式で割り算を行うことは我々が解こうとしている問題と同じくらいの演算コストとなってしまいます．

The idea of using the scale is that it avoids doing any extra computations at the full current precision.  For example, doing the division in the derived inequality above would be as expensive as the problem we're trying to solve!

scaleを用いることは近似に伴う切捨てが許容されるときのみ適当であるということにご注意ください．厳密に条件を満たす必要がある場合，最大精度でその条件を確認する必要があります．

Note that it is important to keep in mind, that using scale is only appropriate where an approximate cutoff is acceptable.  If you need a condition to satisfied strictly, you should check that condition at full precision.

これらが正しく動作するかどうか見てみましょう．

Now lets try it to see if it works

Timing[
      rip  =  AdaptiveNewton[DivisionStep[pin],  N[1/Pi],  
          PrecisionGoal -> 100000]; c]

これは更に改善できます．ここで用いた最終的な精度は終了時点で要求される精度を超えていることにご注意ください．これは精度が高すぎるということを意味しています．乗算コストは演算精度とともに急増するので，これは全体の所要時間に影響を及ぼします．これを解決するためには要求される精度で終了できるように注意深く精度の順番を選ぶ必要があります．

It turns out we can do better yet.  Notice that the final precision used exceeds the goal precision.  This means we have used too much precision.  Since multiplication cost increases quickly with precision, this will affect the overall timing.  A solution is to choose the sequence of precision carefully so we wind up with exactly the goal precision.

AdaptiveNewton[step_,  start_,  
        opts___]  :=  
Module[{goal,  
          prec}, 
c  =  0; 
        goal  =  
          PrecisionGoal  /.  
            Flatten[{opts,  
                Options[AdaptiveNewton]}]; 
prec  =  
          Min[Precision[start],  16]; 
prec  =  
          2^(Floor[Log[2. ,  prec]]  +  
                FractionalPart[
                  Log[2. ,  goal]]); 
SetPrecision[
          NestWhile[AdaptiveNewtonStep[step,  prec],  start,  
            Function[prec  <  2*goal]],  
          goal]
]

これにより評価の1過程を省略することができます．これは最も演算コストがかかるものなので，相当の時間が節約できます．

Doing this saves one extra evaluation, and since it was the most expensive one, a substantial proportion of the time is saved.

Timing[
      rip  =  AdaptiveNewton[DivisionStep[pin],  N[1/Pi],  
          PrecisionGoal -> 100000]; c]

これはMathematica の内部割り算アルゴリズムと比較して優れています．

This compares favorably with Mathematica's internal division algorithm

Timing[1/pin;]

この精度では乗算コストは少し増加するだけです．

It is worth noting that this is only a small multiple of the cost of a multiplication at this precision:

Timing[rip*rip;]

これは単純な計算ですが，精度を少し調整するだけでパフォーマンスが大きく改善されることが分かります．より複雑な計算においては，precision及びaccuracyの調整はあまり明確でないこともありますが，これも同様に重要な場合もあります．

While this is quite a simple calculation, you have seen that some extra thought about controlling precision can bring dramatically better performance.   In more complicated calculations, the control of precision and accuracy is often much less clear, but can be just as important.

また，もう1つの重要な点として，計算速度を劣化させることなく一般的なコードを書くことも可能だということが挙げられます．上記のコードはMathematica の割り算コードの速度とほぼ同じです．しかしこれは割り算に限ったことではありません．例えば，SquareRootsを求めるためにこれを用いることができます．

Another important conclusion is that it is also possible to write general code without sacrificing speed.  The code above is almost as fast as Mathematica's internal division code, but it is not special to division.  For example, we could use it to find SquareRoots.

SquareRootStep[a_][
        x_]  :=  (x^2  +  a)/(2 x)

Timing[rip  =  
          AdaptiveNewton[SquareRootStep[2], Sqrt[2. ],  
            PrecisionGoal -> 100000];]

しかし，これは内部アルゴリズムと比較すると極めて遅くなります．

but it is quite slow compared to the internal algorithm

Timing[N[Sqrt[2],  100000];]

これは各反復演算が割り算を含んでいるため，演算コストが高くなります．以下の関数を用いると

This, was expensive because each iteration included a division.  Using the function

a - (1/(x^2))

2乗根の逆数を求めるために，割り算を行わない方法を構築することも可能です．（2による割り算はMathematica においては非常に簡単です．なぜならそれはバイナリ表現において1ビットシフトに相当するからです．）

for the reciprocal of the square root, it is possible to construct a divisionless step (division by 2 is simple for Mathematica because with our binary representation, it is just a bit shift.

ReciprocalSquareRootStep[a_][x_]  :=  
      x*(3  -  a x^2)/2

割り算を使わない方法を用いると，内部アルゴリズムと比較してそれほど遅くはなりません（内部アルゴリズムではある特別な最終ステップを用いることにより，最後の段階で割り算を使いません）．

Using the divisionless step, it is not a lot slower than the internal algorithm (the internal one avoids the division at the end by taking a special last step)

Timing[rip  =  
          AdaptiveNewton[ReciprocalSquareRootStep[2], Sqrt[2. ],  
            PrecisionGoal -> 100000];]

{0.28000000000000114` Second, Null}

その他の関連例：累積積分
More involved example: Cumulative integrals.

以下によく聞かれる質問を示します．関数fが与えられているとき，以下の式

This is a question which comes up from time to time.  If you have a given function, f, how do you compute

(∫_a^x) f(
            s) \[DifferentialD]s ==  c

を満足するxをどのように求めたらよいのでしょうか．

numerically a value of x so that the equation is satisfied.

ここで以下のような簡単な例を考えてみます．

Lets take a simple example

f[x_]  :=  (1/(Sqrt[π])) Exp[-x^2]

これは以下の累積積分を含んでいます．

with cumulative integral

Cf[x_,  opts___]  :=  
      NIntegrate[f[s], {s, -Infinity,  x},  opts]

ここで，実軸上で1を最大値として99.9%を示すxの極限値を考えます．

and look for what x limit gives 99.9% of the full value of 1 over the real line:

q  =  999/1000;

これを以下の図に示します．

graphically, this looks like

Plot[
        Evaluate[{Integrate[f[s], {s, -Infinity, x}], q}], {x, -3, 
          3}];

この計算を行うのに最も簡単な方法は，FindRootを用いることです．

The simplest way to do it is with FindRoot:

root  =  
      x  /. FindRoot[Cf[x]  -  q, {x, 0},  
          Jacobian -> f[x]]

Cf[root]  -  q

ここでxはp桁のprecision（精度）が望ましいとします．

Now suppose we want to know x to p digits of precision.

p  =  10;

FindRootの中のAccuracyGoalオプションをAccuracyGoal->pに変更するだけでよいはずです．

Presumably, all we have to do is change the AccuracyGoal option in FindRoot to AccuracyGoal->p

root  =  
      x  /.  FindRoot[
          NIntegrate[f[s], {s,  -Infinity,  x}]  -  q, {x, 0},  
          Jacobian -> f[x],  AccuracyGoal -> p  +  2]

しかし結果は要求される誤差範囲に納まっていません．

However, a plot shows that our solution does not lie in the requested error interval.

Plot[
        Cf[(1  +  x) root]  -  q, {x, -10^-p, 
          10^-p}];

経験的に，AccuracyGoalをAccuracyGoal->p + 2に増加させることによってこの問題を解決できることが分かります．

Experimentally, you can find that increasing the AccuracyGoal to AccuracyGoal->p + 2 takes care of this problem:

root  =  
        x  /.  FindRoot[
            NIntegrate[f[s], {s,  -Infinity,  x}]  -  q, {x, 0},  
            Jacobian -> f[x],  
            AccuracyGoal -> p  +  2];
Plot[
        Cf[(1  +  x) root]  -  q, {x, -10^-p, 
          10^-p}];

実際の解に当てはめてみましょう．（この例は的確な解をもつように選びました．このような解は，アルゴリズムの詳細について分析する際に有効であることがあります．）

Now compare it to the exact solution. (Yes, this example was conveniently chosen to have an exact solution.  Such a solution is often very useful when trying to work out the details of an algorithm.)

ExactCf[x_]  :=  (Erf[x]  +  1)/2

Plot[
        ExactCf[(1  +  x) root]  -  q, {x, -10^-p, 
          10^-p},  PlotRange -> All];

正確な解のために要求される誤差範囲を超えてしまいました．実際，拡大したプロットから，誤差は10倍大きすぎることが分かります．

Ouch!  The plot shows that we are NOT within the requested error interval for the exact solution.  In fact, an expanded plot shows that our error is a factor of 10 too large.

Plot[
        ExactCf[(1  +  x) root] - q, {x, -10^(1 - p), 
          10^(1 - p)}];

このような数値コマンドを用いるときは，確実に「正確な」解から指定された誤差内に納めるように推定誤差を十分に慎重に用いるということが大切です．

The entire idea of a numerical command like this is to use error estimates conservative enough to be fairly sure we are within the prescribed error interval from the exact solution.

今まで行われたアプローチには明らかに問題があります．

Clearly something is wrong with the approach taken so far.  What is it?

このプロットは非常に明確にこの問題の原因を実証します．

This plot demonstrates the source of the problem quite clearly:

Plot[
        Cf[(1  +  x) root]  -  
          ExactCf[(1  +  x) root], {x, -10^(1  -  p),  
          10^(1  -  p)},  PlotRange -> All];

関数FindRootが求めている根の値は実質関数の大きさと同じscaleでランダムに変化しています．このような状況ではFindRootは正しい値を求める方法がありません．

The value of the function FIndRoot is trying to find a root for is varying in essentially random way on the same scale as the size of the function.  There is no way FindRoot can be expected to find the correct value under these circumstances.

NIntegrateを用いてよりよい結果を求めようとすると，ローカルな変化のscaleを減少させることができます．

If you ask NIntegrate to do a better job, then we can reduce the scale of local variation:

Plot[
        Cf[(1  +  x) root,  AccuracyGoal -> 12,  
            PrecisionGoal -> 12]  -  
          ExactCf[(1  +  x) root], {x, -10^(1  -  p),  
          10^(1  -  p)},  PlotRange -> All];

これは希望の精度で根を求めるのに十分です．

This should be sufficient to get the root to the desired precision

root  =  
      x  /.  FindRoot[
          Cf[x,  AccuracyGoal -> 12,  PrecisionGoal -> 12]  -  
            q, {x, 0},  Jacobian -> f[x],  
          AccuracyGoal -> p +  2]

さらに，プロットが表示されます．

and a plot shows that it was.

Plot[
        ExactCf[(1  +  x) root]  -  q, {x, -10^-p, 
          10^-p},  PlotRange -> All];

つまりこの問題に関しては，NIntegrateにおけるAccuracyGoal及びPrecisionGoal，FindRootにおけるAccuracyGoal，及び根を求めるために必要なprecisionとaccuracyの間のつり合いの問題が存在します．NIntegrateのコストはAccuracy/Precision Goalsの増加に伴って指数関数的に増加するため，パラメータを平衡点に極めて近い値に設定することが重要です．

So, for this problem, we have a balancing act between the AccuracyGoal and PrecisionGoal for NIntegrate, the AccuracyGoal for FindRoot, and the precision or accuracy we need for the root.  Because the cost of NIntegrate goes up exponentially fast with increase Accuracy/Precision Goals, it is important to get parameters quite close to the balance point.

正確な結果が得られたとき，適切な解を求めるのに必要な制約条件を定めるためにグラフを用いることはそれほど難しくはありません．比較するための正確な解がないまま自動的に制約条件を決定するのはより難しい問題です．以下にパラメータの均衡を取り扱う一般的な方法（完全に解決できるわけではないのですが）を示します．

When you have an exact result, it is not terribly difficult to use graphs to determine the constraints needed to obtain an adequate solution.  It is a more challenging problem to determine the constraints automatically with no exact solution to compare with.  Below, I develop a general (though not completely bulletproof) way of handling the balance of parameters.

根近辺での局所感度は傾きで与えられます．\(TraditionalForm\`\(x\_*\)\)の累乗根近傍の局所的線形モデルはf(x) ≅ (x - \(TraditionalForm\`\(x\_*\)\)) \(TraditionalForm\`f' \((\(x\_*\))\)\)であり，xの誤差はおよそx - \(TraditionalForm\`\(x\_*\)\) ≅ \(TraditionalForm\`\(f(x)\)/
            f' \((\(x\_*\))\)\  ≅ \ \(f(x)\)/f' \((x)\)\)になります．もちろん最後の比は明らかに単純なニュートン法です．

The local sensitivity around the root is given by the slope.  The local linear model around the root \(TraditionalForm\`\(x\_*\)\) is f(x) ≅ (x - \(TraditionalForm\`\(x\_*\)\)) \(TraditionalForm\`f' \((\(x\_*\))\)\), and so the error in x is approximately x - \(TraditionalForm\`\(x\_*\)\) ≅ \(TraditionalForm\`\(f(x)\)/
            f' \((\(x\_*\))\)\  ≅ \ \(f(x)\)/f' \((x)\)\).  The last ratio should, of course, look very familiar since it is plain and simple the Newton step!

ご自分でニュートン反復演算ループを書かれるのでしたら，割り算の例で示したように停止条件として適当な境界を用いることができます．現在のところ，FindRootの停止条件としては関数の大きさのみを用いていますので，傾きに応じてAccuracyGoalを調整する必要があります．通常根の妥当な近似値は，その根においてオプションを調整するのに十分な精度の傾きの近似値を与えます．これは2段階の操作で行うことができます．

If you write your own Newton iteration loop, you could use an appropriate bound on this as a stopping condition (much as we did for the division example).  Currently, since FindRoot only uses the size of the function as a stopping condition, you need to adjust the AccuracyGoal according to the slope.  Typically a reasonable approximation to the root will give you a good enough approximation to the slope at the root to adjust the option.  This can be done with a two step operation.

crude  =  
        x  /.  FindRoot[Cf[x]  -  q,  {x, 0},  
            Jacobian -> f[x]];
slope  =  f[crude]

ここで，関数値の傾きに対する比が目標としている精度より低いことが必要になります．特に

Given that, we want the ratio of the function value to the slope to be less than the precision we are aiming for.  In particular

((10^(-FindRootAG))/slope) ≤ 
      x (10^(-PG))

最小条件は等式であり，両辺の対数をとることにより以下の式が得られます．

The minimal condition is at equality and taking logarithms on both sides, we have

-FindRootAG - (Subscript[log, 10])(slope) == (Subscript[log, 10])(
          x) - PG

さらにScaleの定義より，以下が得られます．

Recalling the definition of Scale, you get

FindRootAG  =  p  -  Scale[slope]  -  Scale[crude]

FindRootを再び使用する前に，NIntegrateで何をするべきかを考える必要があります．これは実際は比較的簡単です．積分値と定数値の差のscaleが-FindRootAGになるまで繰り返すようにFindRootを設定すればよいのです．定数は十分な精度をもち（これが満たされていないと，どのような場合でも結果は劣化します），またそれはNIntegrateで要求される最低の精度であると仮定します．実際，根の近似値付近でFindRootに対して関数が滑らかであるためには（また，上記のようにランダムなプロットでないためには），多少手を加えた方がよいかもしれません．もちろんこれは，NIntegrateが誤差制御も見積りに基づいて行っているため，安全策をとるのにも有効です．

Before calling FindRoot again, it is appropriate to think about what we have to require from NIntegrate.  This is actually relatively simple:  we are asking FindRoot to iterate until the difference of the value of the integral and  a constant be of scale equal to -FindRootAG.  Assuming the constant has sufficient precision (if it doesn't you're sunk in any case!), then that is the minimal accuracy which is required from NIntegrate.  Actually, to be sure that the function appears smooth (and not random like the plots above) to FindRoot near the approximate root, it is a good idea to give a little extra.  This is also helpful to put us on the safe side since, of course, NIntegrate is basing its error control on estimates as well.

NIntegrateExtraAccuracy  =  1;

これらすべてを用いて，根の最終結果を以下のように計算することができます．

With all this, we can compute our final root with

root  =  
      x  /.  FindRoot[
          Cf[x,  AccuracyGoal -> 
                FindRootAG  +  NIntegrateExtraAccuracy,  
              PrecisionGoal -> Infinity]  -  q,  {x,  crude},  
          Jacobian -> f[x],  
          AccuracyGoal -> FindRootAG]

また，以下のプロットは正しい解を表示しています．

And from the plot, it looks like we got it!

Plot[
        ExactCf[(1  +  x) root]  -  q, {x, -10^-p, 
          10^-p},  PlotRange -> All];

これですべての難題が解決しました．これらの考えはすべて統合して関数にする価値があります．そこで，コマンドを融合するパッケージを構築しました．

Now that we have done all of this head scratching, it's worth putting together a function which combines these ideas.  I have put together a package which merges the commands,

FindIntegralLimit1

Needs["FindIntegralLimit1`"]

?FindIntegralLimit1

このパッケージで用いる議論すべきユーティリティは，パッケージAutomaticPGAG.mに要約されています．このパッケージはMathematica の組込み数値関数が用いるオプションAccuracyGoal->AutomaticとPrecisionGoal->Automaticのデフォルトを模倣しています．

A utility this package uses which is worth discussing is encapsulated in the package AutomaticPGAG.m, which mimics the defaults for options AccuracyGoal->Automatic and PrecisionGoal->Automatic which Mathematica's built-in numerical functions use.

単純な関数は操作の中で用いるのが簡単だという点で優れています．例えば以下のようなプロットでは，

An advantage of a nice simple function is that it is easy to use it in operations, such as plotting:

Plot[
        x  /.  FindIntegralLimit1[f[s], {s, -Infinity,  x,  0},  
            q],  {q, 1/2, 9999/10000}];

期待通りの結果を示します．さらに興味深いのは，誤差のプロットです．

Shows the expected behavior.  What is more interesting is a plot of the error:

Plot[
        Log[10, Abs[
            InverseErf[0,  
                2 q  - 1]  -  (x  /.  
                  FindIntegralLimit1[f[s], {s, -Infinity,  x,  0},  
                    q])]],  {q, 1/2, 9999/10000},  
        PlotRange -> All];

プロットの構成は誤差制御の処理を示すものです．これは直接的にはこれから解く問題の特徴ではありません．

The structure of the plot is an artifact of the way error control is handled.  It is not directly a feature of the problem being solved.

これを作るのに多少の時間を要します．精度が要求されるときにコマンドのスピードはより大きな問題となります．

You may have noticed that this took a while to produce.  The speed of the command becomes more of an issue when more precision is requested.

Timing[
      FindIntegralLimit1[f[s], {s, -Infinity,  x,  0},  999/1000,  
        WorkingPrecision -> 24]]

スピードを上げるための1つの方法は，NIntegrateの作業が減るように以下のように入力を書き直すことです．

One way this can be made faster is to rewrite the input so that NIntegrate has less work to do:

Timing[
      FindIntegralLimit1[f[s], {s, Infinity,  x, 2.185124219133004},  
        999/1000  -  1,  WorkingPrecision -> 24]]

しかしこれはいつでも可能なわけではありません．次のセクションではこれを行ういくつかの可能性を扱います．

but that is not always possible.  The next section addresses some possibilities for doing this.

スピードアップの可能性
Some possible speedups

キャッシュによる積分計算値
Cache computed integral values

積分範囲が狭いと，NIntegrateがもっと速くなることは明らかです．NIntegrateが計算を行う領域を限定すれば，時間を節約することができます．

It is clear that doing the smaller integral was much faster for NIntegrate.  By limiting the extent of the region over which NIntegrate has to work, we should be able to save some time.  Because

(∫_a^x) f(
            s) \[DifferentialD]s == (∫_a^((Subscript[x, 1])))f(
              s) \[DifferentialD]s + (∫_((Subscript[x, 1]))^x) f(
              s) \[DifferentialD]s

が成り立つので，以下のように行えます．

We can do the following

cachex  =  -Infinity;
cachei  =  0;
Cfcache[x_,  
          opts___]  :=  (cachei  +=  
            NIntegrate[f[s], {s, cachex,  x},  opts];  cachex  =  x;  
          cachei);

Timing[Cf[0,  WorkingPrecision -> 24]]

Timing[Cfcache[0,  WorkingPrecision -> 24]]

Timing[
      Cf[2.185124219133004264123974081742574`24,  
        WorkingPrecision -> 24]]

Timing[
      Cfcache[2.185124219133004264123974081742574`24,  
        WorkingPrecision -> 24]]

しかし，これはあまり正確でないので，この方法は注意深く行わなければなりません．NIntegrateはグローバル誤差の見積りを計算します．NIntegrateでいくつかの結果を追加すると，これは正確ではなくなります．しかし，付け加えられた数があまり多くなければ，付加的な誤差はあまり大きくはなりません．例えば，10個のインターバルを統合した場合，NIntegrateExtraPrecisionに1個付加した場合と同様の精度の結果が期待できます．

You have to be careful with this approach, however because it is not quite as accurate.  NIntegrate computes a global error estimate.  This is lost  when you add the results of several invocations in NIntegrate.   However, if the number added together is not too large, the extra error is not too large.  For example, even if you add 10 intervals together, we can expect the same quality result by adding 1 extra to NIntegrateExtraPrecision.

これをこのパッケージの第2バージョンFindIntegralLimit2.mに付け加えました．

I've put this into a second version of the package, FindIntegralLimit2.m.

Needs["FindIntegralLimit2`"]

これはこの例題でより速く動作します．

and, sure enough, it is faster on this example.

Timing[
      FindIntegralLimit2[f[s], {s, -Infinity,  x,  0},  999/1000,  
        WorkingPrecision -> 24]]

NDSolveを用いる
Use NDSolve

NDSolveもまた数値積分に用いることができますが，グローバル誤差制御を備えていないため通常は適切ではありません．しかし，大雑把な見積りを行う場合は非常に適切に動作することがあります．

NDSolve can be used as a numerical integrator also, though since it does not have global error control, it is usually not appropriate.  For the crude approximation, however, it can in some cases work quite nicely.

NDSolveは無限の極限値を扱うことができないため，以下の2つのどちらかを行わなければなりません．1つはNIntegrateが自動的に行うように変数変換を行うということです．これはそれほど難しくはありませんが，ここでは取り扱いません．

Since NDSolve cannot handle infinite limits, there are two things we can do.  One is to do a variable transformation much like what NIntegrate does automatically.  While this is not exceptionally difficult, I will not cover this method here.

最も簡単なのは，計算に用いる積分の有限部分を特定することです．NDSolveがどのように動作するかを見るために，最小値が0になるように問題を変えてみましょう．

The easiest is to identify a finite part of the integral to use.  So that you can see how NDSolve works, lets change the problem so the lower limit is 0.

ndsol  =  
        First[u  /.  
            NDSolve[{u'[s]  ==  f[s],  u[0]  ==  0},  
              u,  {s, 0,  10}]];

Plot[ndsol[s], {s, 0, 10},  PlotRange -> All];

FindRoot[ndsol[s]  -  499/1000, {s, 1, 0, 10}]

これをより効率的に行うために，NDSolveはStoppingTestと呼ばれるオプションを備えています．これは指定の条件がTrueとなったときにNDSolveを停止させるものです．

To make this a little more efficient, NDSolve has a neat little option called StoppingTest, which causes NDSolve to stop when a given condition evaluates to True.

ndsol  =  
      First[u  /.  
          NDSolve[{u'[s]  ==  f[s],  u[0]  ==  0},  
            u,  {s, 0,  10},  
            StoppingTest -> (u[s]  >  499/1000)]]

これは最後のステップ s == 2.32251で，uが閾値に達したことを表しています．したがって理想的には，FindRootの検索範囲を

This means that in the last step to s == 2.32251, u crossed the threshold.  So ideally, we want to narrow the FindRoot search to between

ndend  =  ndsol[[3, 1, -1]]

及び

and

ndprev  =  ndsol[[3, 1, -2]]

の範囲に狭めた方がよいと考えられます．FindRootが呼び出されたとき，異なるブレント（Brent）のアルゴリズムが用いられます．これにより根はカッコの中に表示されます．

When FindRoot is called, it will use a variant of Brent's algorithm which keeps the root bracketed.

FindRoot[ndsol[s]  -  499/1000,  {s,  ndprev,  ndend}]

グローバル誤差制御を備えていないということから，NDSolveを最終段階の近似には用いたくありません．しかし，これは大雑把に近似を行うには非常に適しています．NDSolveの使用法はパッケージの第3バージョン FindIntegralLimit3.mに組み込まれています．

Because of the fact that NDSolve does not have global error control, we do not want to use it for the final approximation.  However, this is quite adequate for the crude approximation.  The use of NDSolve has been incorporated into a third version of the package, FindIntegralLimit3.m

Needs["FindIntegralLimit3`"]

Timing[
      FindIntegralLimit3[f[s], {s, -Infinity,  x,  0, 10},  
        999/1000,  WorkingPrecision -> 24]]

初期の根の分離がより難しい場合にこの方法を用いると，より時間を節約できます．

This approach could save more time for an example where the initial root isolation was more difficult.

FindIntegralLimit

以下に今までのバージョンに含まれていないもので最終パッケージが備えているべきものを示します．

So what should the final package have that I have not included in any of my three versions so far.

• 大雑把な近似値を求める方法を切り換えるきちんとした方法が必要です．1つの可能性として提案されているのは，FindIntegralLimit3で必要な引数を追加するということです．これの欠点としては単純に引数を追加することは直観的にはNDSolveを使うように変更する明らかな方法ではないということです．

• It should have a tidy way of switching methods for finding the crude approximation.  One possibility is suggested by the extra argument required by FindIntegralLimit3.  A drawback of this is that simply adding an extra argument is not intuitively an obvious way to change to using NDSolve.

• 初期条件の精度が機械精度より高い場合は，FindRootの動作と同様に，大雑把な近似は初期条件の精度で行われるべきです．これにより，機械精度の数値に不要情報を与えてしまう関数をお持ちなら，意味のない計算に時間を浪費しなくてすみます．

• When the precision of the starting condition is bigger than machine precision, the crude approximation should likely be done with that precision, similar to the way FindRoot works.  That way, if you have a function which gives garbage for machine numbers, you are not wasting time with useless calculations.

• Nintegrate，NDSolve，及びFindRootの誤動作を注意深く確認し，趣味のよいメッセージを表示して優雅に終了することが必要です．

• It should check carefully for failures in NIntegrate, NDSolve  and FindRoot and exit gracefully with a tasteful message.

• （FindRootの）計算の最終段階で何度か精度を高める必要がある場合，精度が変更になるたびに自動的に積分計算をやり直す方法があるとよいでしょう．（バージョン2と3では，FindRootが精度を高めてしまう問題を回避するために，FindRootを最後に呼び出すときに最高の精度を用います．NIntegrateが使用できる実用的な精度を限定するため，このような精度の大きな増加はまれです．）

• When the precision needs to be increased several times in the last step (done with FindRoot), it would be better to have an automatic way to restart the computation of the integral from scratch whenever the precision is changed.  (In versions 2, 3 we go to full precision for the last call to FindRoot to avoid the problem of FindRoot increasing precision.  Such a large increase in precision will be rare since NIntegrate limits the practical precision which can be used.

• 大雑把な近似において（ヤコビの）傾きがゼロ付近であるとき，別の方法を用いた方がよいでしょう．

• It would be good to come up with alternative schemes in the case where the slope (Jacobian) is zero or nearly so at the crude approximation.

• テストを行いましょう．ここでは比較的簡単な例を用いて関数の考え方を示しました．確実なパッケージを作り上げるには，種々のテストを行わなければならないでしょう．

• Tests.  I have just shown the idea of the function on a relatively simple example.  To get a robust package working, you would have to come up with a variety of tests.

最後に
Summary

accuracyとprecisionの概念をよく理解したら，AccuracyGoalとPrecisionGoalオプションの値を特定したときに，数値関数からいかによく解を予想できるかがお分かりになるでしょう．またこれにより，さらにいくつかの基本的な数値操作を結合して，より複雑な数値関数を書くことができます．

Once you have a good understanding of the concepts of accuracy and precision, you can gain a good understanding of how good you should expect solutions from numerical functions to be when you specify values for the AccuracyGoal or PrecisionGoal options.  From this, it is possible to write more complex numerical functions which combine several mroe elementary numerical operations.