これは、問題のもっともらしい方法について尋ねた相互検証済みの質問への続きです。この質問はよりプログラミング指向であるため、SO に投稿します。
バックグラウンド
既知の日付が 1 年以上にわたる曲線があります。この曲線の y 値は、毎日の気温と塩分の記録から計算されたd18O 値の予測です。また、炭酸カルシウムからなるシェルの d18O 値も測定しました。これらの値は、距離軸に沿って測定されます。最初と最後の測定は、曲線の開始と終了とほぼ同時に (正確ではありませんが) 行われます。
d18O 値は、未知の確率誤差内で曲線の予測値と一致することが知られています。測定値の x 軸を変更することによって (または少なくともインデックスを曲線のインデックスと一致させることによって)、測定値を曲線に最適に適合させたいと考えています。このようにして、測定値の日付を推定し、年間の殻の成長率をさらに推定することができます。成長率は変動すると予想され、成長の停滞（成長の停止）が生じる可能性があります。ただし、測定値間の増加は > 0 でなければなりません(制約)。
サンプルデータ
データセットの例 (curveおよびmeasured) は次のとおりです。
meas <- structure(list(index = 1:10, distance = c(0.1, 1, 3, 5, 7, 8,
13, 20, 22, 25), value = c(3.5, 4.2, 4.5, 4.4, 4.7, 4.8, 5.1,
4.9, 4.1, 3.7)), .Names = c("index", "distance", "value"), class = "data.frame",
row.names = c(NA, -10L))
curve <- structure(list(date = structure(c(15218, 15219, 15220, 15221,
15222, 15223, 15224, 15225, 15226, 15227, 15228, 15229, 15230,
15231, 15232, 15233, 15234, 15235, 15236, 15237, 15238, 15239,
15240, 15241, 15242, 15243, 15244, 15245, 15246, 15247, 15248,
15249, 15250, 15251, 15252, 15253, 15254, 15255, 15256, 15257,
15258, 15259, 15260, 15261, 15262, 15263, 15264, 15265, 15266,
15267, 15268, 15269, 15270, 15271, 15272, 15273, 15274, 15275,
15276, 15277, 15278, 15279, 15280, 15281, 15282, 15283, 15284,
15285, 15286, 15287, 15288, 15289, 15290, 15291, 15292, 15293,
15294, 15295, 15296, 15297, 15298, 15299, 15300, 15301, 15302,
15303, 15304, 15305, 15306, 15307, 15308, 15309, 15310, 15311,
15312, 15313, 15314, 15315, 15316, 15317, 15318, 15319, 15320,
15321, 15322, 15323, 15324, 15325, 15326, 15327, 15328, 15329,
15330, 15331, 15332, 15333, 15334, 15335, 15336, 15337, 15338,
15339, 15340, 15341, 15342, 15343, 15344, 15345, 15346, 15347,
15348, 15349, 15350, 15351, 15352, 15353, 15354, 15355, 15356,
15357, 15358, 15359, 15360, 15361, 15362, 15363, 15364, 15365,
15366, 15367, 15368, 15369, 15370, 15371, 15372, 15373, 15374,
15375, 15376, 15377, 15378, 15379, 15380, 15381, 15382, 15383,
15384, 15385, 15386, 15387, 15388, 15389, 15390, 15391, 15392,
15393, 15394, 15395, 15396, 15397, 15398, 15399, 15400, 15401,
15402, 15403, 15404, 15405, 15406, 15407, 15408, 15409, 15410,
15411, 15412, 15413, 15414, 15415, 15416, 15417, 15418, 15419,
15420, 15421, 15422, 15423, 15424, 15425, 15426, 15427, 15428,
15429, 15430, 15431, 15432, 15433, 15434, 15435, 15436, 15437,
15438, 15439, 15440, 15441, 15442, 15443, 15444, 15445, 15446,
15447, 15448, 15449, 15450, 15451, 15452, 15453, 15454, 15455,
15456, 15457, 15458, 15459, 15460, 15461, 15462, 15463, 15464,
15465, 15466, 15467, 15468, 15469, 15470, 15471, 15472, 15473,
15474, 15475, 15476, 15477, 15478, 15479, 15480, 15481, 15482,
15483, 15484, 15485, 15486, 15487, 15488, 15489, 15490, 15491,
15492, 15493, 15494, 15495, 15496, 15497, 15498, 15499, 15500,
15501, 15502, 15503, 15504, 15505, 15506, 15507, 15508, 15509,
15510, 15511, 15512, 15513, 15514, 15515, 15516, 15517, 15518,
15519, 15520, 15521, 15522, 15523, 15524, 15525, 15526, 15527,
15528, 15529, 15530, 15531, 15532, 15533, 15534, 15535, 15536,
15537, 15538, 15539, 15540, 15541, 15542, 15543, 15544, 15545,
15546, 15547, 15548, 15549, 15550, 15551, 15552, 15553, 15554,
15555, 15556, 15557, 15558, 15559, 15560, 15561, 15562, 15563,
15564, 15565, 15566, 15567, 15568, 15569, 15570, 15571, 15572,
15573, 15574, 15575, 15576, 15577, 15578, 15579, 15580, 15581,
15582, 15583, 15584), class = "Date"), index = 1:367, value = c(3.33,
3.35, 3.36, 3.38, 3.4, 3.42, 3.43, 3.45, 3.47, 3.48, 3.5, 3.52,
3.53, 3.55, 3.56, 3.58, 3.6, 3.61, 3.63, 3.64, 3.66, 3.67, 3.69,
3.7, 3.72, 3.73, 3.75, 3.76, 3.78, 3.79, 3.81, 3.82, 3.83, 3.85,
3.86, 3.88, 3.89, 3.9, 3.92, 3.93, 3.94, 3.96, 3.97, 3.98, 3.99,
4.01, 4.02, 4.03, 4.04, 4.06, 4.07, 4.08, 4.09, 4.1, 4.11, 4.13,
4.14, 4.15, 4.16, 4.17, 4.18, 4.19, 4.2, 4.21, 4.22, 4.23, 4.24,
4.25, 4.26, 4.27, 4.28, 4.28, 4.29, 4.3, 4.31, 4.32, 4.33, 4.33,
4.34, 4.35, 4.36, 4.36, 4.37, 4.38, 4.38, 4.39, 4.4, 4.41, 4.41,
4.42, 4.42, 4.43, 4.44, 4.44, 4.45, 4.45, 4.46, 4.46, 4.47, 4.47,
4.47, 4.48, 4.48, 4.49, 4.49, 4.49, 4.5, 4.5, 4.5, 4.51, 4.51,
4.51, 4.52, 4.52, 4.53, 4.53, 4.53, 4.54, 4.54, 4.54, 4.55, 4.55,
4.56, 4.57, 4.57, 4.58, 4.58, 4.59, 4.6, 4.61, 4.61, 4.62, 4.63,
4.64, 4.64, 4.65, 4.66, 4.67, 4.67, 4.68, 4.69, 4.7, 4.7, 4.71,
4.72, 4.72, 4.73, 4.74, 4.74, 4.75, 4.75, 4.75, 4.76, 4.76, 4.76,
4.76, 4.76, 4.76, 4.76, 4.76, 4.76, 4.75, 4.75, 4.75, 4.75, 4.74,
4.74, 4.73, 4.73, 4.73, 4.72, 4.72, 4.72, 4.71, 4.71, 4.71, 4.71,
4.7, 4.7, 4.7, 4.71, 4.71, 4.71, 4.71, 4.72, 4.72, 4.73, 4.74,
4.75, 4.75, 4.76, 4.78, 4.79, 4.8, 4.81, 4.82, 4.83, 4.84, 4.85,
4.86, 4.88, 4.89, 4.9, 4.91, 4.92, 4.92, 4.93, 4.94, 4.95, 4.95,
4.95, 4.96, 4.96, 4.96, 4.96, 4.96, 4.95, 4.95, 4.95, 4.94, 4.93,
4.92, 4.92, 4.91, 4.9, 4.89, 4.88, 4.87, 4.86, 4.85, 4.84, 4.83,
4.82, 4.8, 4.79, 4.78, 4.77, 4.76, 4.75, 4.75, 4.74, 4.73, 4.72,
4.72, 4.71, 4.71, 4.71, 4.7, 4.7, 4.7, 4.7, 4.7, 4.7, 4.7, 4.7,
4.7, 4.7, 4.7, 4.7, 4.7, 4.69, 4.69, 4.69, 4.69, 4.69, 4.69,
4.69, 4.69, 4.68, 4.68, 4.68, 4.67, 4.67, 4.67, 4.66, 4.65, 4.65,
4.64, 4.63, 4.62, 4.61, 4.6, 4.59, 4.58, 4.57, 4.56, 4.55, 4.54,
4.53, 4.51, 4.5, 4.49, 4.48, 4.47, 4.46, 4.45, 4.43, 4.42, 4.41,
4.4, 4.39, 4.38, 4.37, 4.36, 4.35, 4.34, 4.33, 4.32, 4.32, 4.31,
4.3, 4.29, 4.28, 4.28, 4.27, 4.26, 4.25, 4.24, 4.24, 4.23, 4.22,
4.21, 4.21, 4.2, 4.19, 4.18, 4.17, 4.17, 4.16, 4.15, 4.14, 4.14,
4.13, 4.12, 4.12, 4.11, 4.1, 4.09, 4.08, 4.08, 4.07, 4.06, 4.05,
4.05, 4.04, 4.03, 4.02, 4.02, 4.01, 4, 4, 3.99, 3.98, 3.97, 3.97,
3.96, 3.95, 3.94, 3.94, 3.93, 3.92, 3.92, 3.91, 3.9, 3.9, 3.89,
3.88)), .Names = c("date", "index", "value"), row.names = c(NA,
-367L), class = "data.frame")
...そして、これがどのように見えるかです：
library(ggplot2)
library(scales)
library(gridExtra)
p.curve <- ggplot() + geom_line(data = curve, aes(x = date, y = value)) + scale_x_date(name = "Month", breaks = date_breaks("months"), labels = date_format("%b")) + labs(title = "curve")
p.meas <- ggplot(meas, aes(x = distance, y = value)) + geom_point(color = "red") + labs(title = "measured", x = "Distance (mm)")
grid.arrange(p.curve, p.meas, ncol = 1)
実際の問題
の x 軸を変更することにより、measRが適合する数学的/統計的方法を見つけたいと考えています。さらに、適合した「x軸」を相互に比較するために、ある種の適合度統計を取得したいと考えています（異なる制約で複数のモデルを実行する場合）。私は「x 軸モデル」を成長モデルと呼んでいます。値間の距離が > 0 でなければならないことを指定して、フィッティングを制限したい。つまり、 の値は の値の後に発生する必要がある。また、成長率 (つまり、隣接する 2 つのインデックス ポイント間の最大距離) を制限できるようにしたいと考えています。これを実証するために、手動で行います。curvemeasmeasMeasindex == 2index == 1
ggplot() + geom_line(data = curve, aes(x = index, y = value)) + geom_line(data = meas, aes(x = index, y = value), color = "red", linetype = 2) + scale_x_continuous(breaks = seq(0,370,10)) + scale_y_continuous(breaks = seq(3,5,0.1))
まず、meas(赤い破線)のインデックスのいくつかをcurve(黒い線) のインデックスに固定する必要があります。最初と最後のポイントに加えて、最も高い値を持つポイントを固定することを選択します。
anchor <- data.frame(meas.index = c(1,7,10), curve.index = c(11,215,367))
example.fit <- merge(meas, anchor, by.x = "index", by.y = "meas.index", all = T, sort = F)
example.fit <- example.fit[with(example.fit, order(distance)),]
次に、これらの固定点の間で直線的な成長を仮定します。成長はcurve指数に沿って進みます。Curveには 1 日 1 つの値があります。したがって、成長は毎日の規模になります。
library(zoo)
example.fit$curve.index <- round(na.approx(example.fit$curve.index),0)
この後、インデックスを日付に置き換えて結果をプロットします。
library(plyr)
example.fit$date <- as.Date(mapvalues(example.fit$curve.index, from = curve$index, to = as.character(curve$date)))
a <- ggplot() + geom_line(data = curve, aes(x = date, y = value)) + geom_point(data = example.fit, aes(x = date, y = value), color = "red") + scale_x_date(limits = range(curve$date), name = "Month", breaks = date_breaks("months"), labels = date_format("%b"))
b <- ggplot(example.fit, aes(x = date, y = distance)) + geom_line() + scale_x_date(limits = range(curve$date), name = "Month", breaks = date_breaks("months"), labels = date_format("%b"))
grid.arrange(a,b)
上のプロットは、3 つのアンカー ポイントに基づく結果の適合を示しています。以下のプロットは、1 日間隔で時間に沿ったモデル化された成長を示しています。3 月初めの成長曲線の曲がりは、私には理解できない面白い数学的アーティファクトです (パッケージのna.approx関数による)。zoo
私は何を試しましたか
以前の質問から、ダイナミック タイム ワーピングが解決策になる可能性があることがわかりました。また、dtw 関数を含むR パッケージも見つけました。良い。実際、動的タイムワーピングは、その質問のサンプルデータセットで機能しましたが (制約の設定を除く)、 (と呼ばれる)curveよりもはるかに多くのデータポイントがあるこのデータセットでは機能しません。measpoints前の質問で）。スペースを節約するために、ここにコード/図をコピーしません。その質問に対する私の答えで、私が試したことを見ることができます。問題は、最も単純なパターンを除いて、どのステップ パターンもこれらのタイプのデータを処理できないことです。最も単純なステップ パターンは、測定値を曲線に数回一致させます。これは避けたいことです。測定ポイントごとに定義された日付が必要だからです。また、測定ポイント間で成長率が 0 を超えなければならないという制約を設定するのは難しいようです。
質問
私の質問は 2 つあります。まず、動的タイム ワーピングよりも問題を解決するためのより良い方法はありますか? 第二に、Rで実際にこれを行うにはどうすればよいですか? .
EDITS 9. 2013年12月私は質問をより明確にしようとしました。